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1. Introduction 

At the time of the writing of this review, the tenth anniversary of the discovery of the 
first extrasolar planet around a solar-type star is being celebrated. When this detection 
was announced, it came as a shock, not because an extrasolar planet was at last found 
around a star like our Sun — nobody really had doubts that such a detection would occur 
one day or another — but because of the very short period of this planet, 51 Pegasi. It 
orbits its parent star in only 4 days! It is 100 times closer to the star than Jupiter is 
from the Sun. The temperatures at these distances from the star are so high that the 
material needed to assemble the core of such planets cannot be found in a solid phase 
there. The implication is that the planet formed further out and moved in. It did not 
take long for the theorists to propose an explanation: tidal interactions between the 
planet and the disk in which it formed led to orbital migration. The theory of such 
processes had already been developed, about 15 years before the discovery of 51 Pegasi, 
and orbital migration of planets had been predicted. The real surprise, then, was that 
nobody ever predicted that planets might be found on very short period orbits. 

As of today, about 150 planets have been detected. Their orbital chatacteristics 
suggest that planct/disk tidal interactions is very important in determining the 
characteristics of planetary systems. This review is divided in four sections. The first 
section is devoted to the observations of extrasolar planets. We first present the different 
methods used to detect planets, and review the statistical properties of the systems 
observed so far. In the second section, we review the theories of planet formation. 
We describe the different stages that lead to the build-up of a terrestrial planet or a 
planetary core (grain sedimentation, obtention of planetesimals, runaway accumulation), 
and the process by which an enveloppe of gas may be captured to produce a giant 
planet. We also discuss giant planet formation by gravitational instabilities. The third 
section is devoted to disk/planet interactions and the different types of orbital migration. 
Finally, in the fourth section, we review some results about multiplanet systems and 
their interactions. 

2. Properties of exoplanet systems 

The first object outside our solar system with a mass in the range of planetary masses 
was detected in 1992 around the millisecond pulsar PSR1257+12a (Wolszczan & Frail 
1992). Three planets of 0.02, 4.3 and 3.9 Earth masses orbiting at 0.19, 0.36 and 0.46 
astronomical units (AU) from the pulsar, respectively, were actually found. Another 
object of about 100 Earth masses with a separation of 40 AU was reported more recently 
(Joshi & Rasio 1997, Wolszczan 1996). In 1993, Backer et al announced the detection 
of a 2.5 Jupiter mass (Mj) planet orbiting the miUisecond pulsar 1620-26 at a distance 
of 23 AU. This pulsar is in a binary system with a white dwarf companion (Thorsett et 
al 1999) and a separation of about 40 AU. It is not clear whether these planets formed 
before the explosion of the parent supernova or after, either from the material which 
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was ejected during the explosion or, in the case of the pulsar 1620-26, in the disk of 
material transferred from the white dwarf. 

The first extrasolar planet around a sun-like star (51 Pegasi b) was reported 
in 1995 by Mayor & Queloz. It was detected from Haute-Provcnce Observatory, 
and immediately confirmed by Marcy & Butler (1995), who observed it from Lick 
Observatory (see the review by Marcy & Butler 1998 for more details about the detection 
of the first eight extrasolar planets) . Since then, there has been an uninterrupted flow of 
detections. As of today (16 March 2005), 152 planets in 134 planetary systems (including 
14 multiple planet systems) have been reported around solar type stars. The hghtest 
object has a mass of 0.042 Mj. We will focus here on the objects with masses smaller 
than 13 Mj, as heavier objects are usually considered to be brown dwarfs. This review is 
concerned with extrasolar planets around solar type stars, so we will not discuss further 
planets around pulsars. 

2.1. Detection methods 
2.1.1. Radial velocity: 

So far, almost all the planets have been detected by this technique, which consists in 
measuring the Doppler shift due to the motion of the star around the center of mass of 
the system. The velocity of the star projected along the fine of sight (or radial velocity) 
is, to within a constant, v — Kf{t), where the function f{t) varies periodically with 
time t around the orbit and has zero mean and 



Here G is the constant of gravitation, T is the orbital period, is the mass of the 
star, Mp is the planet mass, i is the angle of the fine of sight with respect to the 
perpendicular to the orbital plane and e is the eccentricity of the orbit. The observation 
of as a function of time gives T directly. If the orbit is circular, f{t) is a sinusoid 
with unit amplitude. For non circular orbits the amplitude can reach 1 + e and the 
departure from a sinusoid allows determination of the eccentricity e. Given T and e, 
one can deduce M^sini from the measurement of K (Mj, being known). Note that since 
Mp -C M*, the semi-major axis a can be calculated from T. 

Since i is unknown, only the lower limit MpSini for the planet mass, or projected 
mass, can be obtained. This was of course an important matter of debate when the very 
flrst planet was detected. However, for a population of about 150 objects, the ratio of 
the real mass to the projected mass is on average on the order of unity. The expression 
for K shows that only relatively massive planets with short periods can be detected 
with the radial velocity method. Using Mp <^ M^, we can rewrite K and T in the form: 
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The detection limit of the instruments that are currently being used is about 
1 m s~^. In principle, Jupiter mass planets could be detected rather far from their 
parent star. However, the detection is confirmed only if the planet if monitored for at 
least one orbital period. Given that surveys began about 10 years ago, only planets 
with semi-major axes less than about 4-5 AU can be observed. Note that a Jupiter 
mass planet with a = 5 AU around a solar mass star induces a velocity f ~ 11 m s~^. 

For an Earth mass planet with a = 1 AU around a solar mass star, w ~ 0.1 m s^^. 
It is not clear whether such a precision could be reached. Such a signal would have to 
be extracted from the "noise" at the stellar surface, i.e. random Doppler shifts due to 
atmospheric oscillations, stellar spots, or long period effects caused by magnetic cycles. 
Note that the distance to the star is not a hmiting factor as long as enough photons can 
be collected. 

Several long-term radial velocity surveys for planets are ongoing in different 
countries (for a list see Ferryman et al. 2005). They all use ground-based telescopes. 
New detections are regularly announced. So far all the planets found with this technique 
are within about 100 parsecs (pc) from the Earth. 

2.1.2. Transits: 

The first transit of an extrasolar planet (HD 209458 b) was detected in September 
1999 (Charbonneau et al 2000). The planet, which is on an orbit with a = 0.045 AU, 
had already been observed by the radial velocity method. The detection of the transit, 
in addition to proving beyond doubt that the Doppler shift was indeed caused by a 
planet, enabled the planet mass (0.63 Mj) to be determined, as sini has to be close 
to 1 for a transit to be seen. It also enabled the planet radius (1.3 Jupiter radius) to be 
measured, as the relative variation of the stellar flux during the transit is: 



where Rp is the planet radius and i?^ is the stellar radius. For Jupiter orbiting the Sun, 
AF^/F^ = 1%, whereas for the Earth it is 0.01%. From the ground, the accuracy is 
about 0.1% and is limited by the variable extinction of the atmosphere. Prom space, 
the accuracy reaches 10~^ % so that terrestrial planets could be detected without any 
problem. However, the probability of detecting a transit, B^/a, is rather small. For 
— 1 Rq, this probability is about 10% for a — 0.05 AU, i.e. for the planets on the 
tightest orbits known to date. Therefore, a large number of stars have to be observed 
for a significant number of transits to be detected. 

Note that the mean mass density of HD 209458 b can be calculated from its mass 
and radius. It is found to be about 400 kg m^'^, significantly smaller than the mean 
density of Saturn (700 kg m~^), the least dense of the planets of our solar system. This 




(4) 



Planet formation and migration 



5 



confirms that this planet is made mainly of gas. The atmosphere of the planet was 
further detected by HST, wliich observed additional sodium absorption due to the light 
passing through the planet atmosphere as it transited across the star (Charbonneau et 
al 2002). More recently, atomic hydrogen, oxygen and carbon have also been detected 
in the extended envelope around the planet (Vidal-Madjar et al 2003, 2004), suggesting 
the atmosphere is evaporating. 

The first planet found by the transit method, TrES-1, was detected by the STARE 
(STellar Astrophysics & Research on Exoplanets) program (Alonso et al 2004), and 
confirmed 8 days later by an amateur astronomer in Belgium. Since then, 5 more planets 
have been found by this method by OGLE (Optical Gravitational Lensing Experiment). 
All these planets have lately been confirmed by radial velocity measurements. Their 
separation from the central star ranges from 0.0225 to 0.047 AU. The planets found by 
OGLE are about 1,500 pc away, much further away than the planets found by the radial 
velocity method. 

Several ground-based surveys for planet transits are ongoing. A few space missions 
dedicated at least in part to the planet search are also being developed. The first 
experiment, which will be launched in 2006, is COROT (for COnvection Rotation and 
planetary Transits), a mini-sateUite developed by the French National Space Agency 
(CNES). Its precision is 0.03%, so that it will detect only planets with a radius larger 
than twice the Earth radius. It will monitor about 60,000 stars during 2.5 years. The 
next experiment is KEPLER, which will be launched in 2007 by NASA. Its precision 
is 10~^ % and it will monitor 10^ stars during 4 years. It is expected to find at least 
50 terrestrial planets. 



2.1.3. Astrometry: 



The radial velocity method described above uses the motion of the star along the 
line of sight. It is also possible to measure the motion of the star perpendicular to the 
line of sight and to derive the characteristics of the planet and its orbit from it. If the 
planet is on a circular orbit, the angular motion of the star is: 

Mpo/IAU 

a (arcsec) = irr—rT, > (5) 

where d is the distance to the star, measured in parsecs. Since d and are known and 
a is given by the periodicity of the star motion, Mp can be derived from the measure 
of a. The motion of the Sun caused by Jupiter and seen from a distance of 10 pc is 
500 //as, whereas that due to the Earth, and seen from the same distance, is 0.3 //as. 

^From the ground, the interferometers provide the best astrometric accuracy. It is 
about 1 mas for the VLTI and KECK in the near IR, and should reach 10 /tas in the 
near future over narrow fields of view. Surveys for planets with these interferometers 
will begin in the very near future. 

For space observations, the accuracy will be on the order of 1 /las. Therefore, since 
the stars that are being observed are at least at a distance of a few parsecs, it will not 
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be possible to detect terrestrial planets in habitable zones with this method, which is 
better suited for the observation of giant planets. Note that astrometric measurements 
are easier for longer period systems. This method and the radial velocity method are 
thus complementary. 

To date, only one space mission, Hipparcos (which was launched by ESA in 1989 
and operated until 1993), was dedicated to astrometric measurements. Although it 
was not meant to detect planets, it did observe stars around which planets were later 
observed. Its measurements were therefore used to constrain the mass of some of the 
planets discovered by the radial velocity method. 

Two space experiments are currently being developed. SIM (Space Interferometry 
Mission) is scheduled for launch by NASA in 2010. Its accuracy will be 1 /las. It will 
not observe a very large number of stars but will focus in particular on the multiple 
systems already detected from the ground to get precise measurement of the orbital 
characteristics. GAIA, which is scheduled for launch by ESA in 2011, will have an 
accuracy of 2 — 10 //as and will observe ~ 10^ stars within a distance of about 200 pc 
from the Sun during 5 years. Simulations indicate that GAIA could detect between 
10,000 and 50,000 giant planets. 

2.1.4. Microlensing: 

This technique relies on a background star to provide a source of light. Then, fore- 
ground stars act as gravitational lenses when they pass in front of the background star 
and cause the source star to suddenly brighten by as much as 1000 times. This typically 
lasts for a few weeks or months. While this is the normal pattern of a microlensing event, 
things are substantially different when the lensing (foreground) star has a smaller com- 
panion. If a planet is positioned close enough to the lensing star, it will indeed perturb 
the light curve by adding two successive spikes of brightness, superimposed upon the 
regular pattern of the microlensing event (Mao & Paczynski 1991, Gould & Loeb 1992, 
Bennett & Rhie 1996). Microlensing is most sensitive to planets at a separation ~ 1- 
5 AU from the lens star. In addition, effects on light curves are large even for planets of 
less than an Earth mass. This technique is therefore well suited for finding terrestrial 
planets in habitable zones. An inconvenient aspect of the microlensing technique is that 
the events do not repeat. 

Several collaborations, from around the world, are presently searching for planets 
with the microlensing technique. Two of them (OGLE and MOA) are survey programs 
which monitor large number of stars. The two other projects (PLANET and MicroFUN) 
pick up on the alerts that are issued by the survey collaborations when a microlensing 
event is detected and do a follow-up monitoring of the event. 

One planet (OGLE235-MOA53) has been found in 2004 simultaneously by the 
MOA and OGLE collaborations (Bond et al 2004). Its mass and semi-major axis are 
about 2 Mj and 3 AU, respectively. It is located at a distance of 5.2 kpc from Earth, 
which is much too far away for the radial velocity method to be used. 
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So far, the microlensing programs have mainly put hmits on the percentage of stars 
harboring planets within a given range of separations. Analysis of 5 years of PLANET 
photometry indicates that less than 33% of the M dwarfs in the Galactic bulge can have 
Jupiter-mass companions within a distance between 1.5 and 4 AU. The upper limit is 
raised to 45% if the range of separations considered is 1-7 AU (Gaudi et al 2002). 

A space-based gravitational microlensing program (GEST, for Galactic Exoplanet 
Survey Telescope) is currently being studied by NASA. GEST, which would be able to 
detect planets with masses as low as that of Mars, would survey more than 100 million 
stars. It would be expected to find about 100 Earth mass planets. 

2.1.5. Direct imaging: 

The direct detections of the transiting planets HD 209458 b (Dcming et al. 2005) 
and TrES-1 (Charbonneau et al 2005) were reported simultaneously in March 2005. 
These were the very first observations of photons emitted by extrasolar planets. The 
detections were made by the Spitzer Space Telescope by observing the decrement in 
infrared flux during the secondary eclipses, when the planets passed behind their parent 
star. These observations enable the brightness temperature of the planets at a given 
wavelength to be calculated and, assuming some model for the atmosphere, a blackbody 
temperature to be estimated. 

In general, direct imaging of extrasolar planets will enable long period objects to 
be detected, and therefore will complement radial velocity surveys. Very low mass non 
transiting objects can currently be detected only if they are isolated (e.g., free floating 
planets). Bound planets are much more difficult to observe because of the very high 
star to planet flux ratio. The contrast is smaller in the IR, which is therefore the most 
favourable range of wavelengths for planet imaging. The light coming from the star can 
in addition be blocked by using a coronograph or nulling techniques on interferometers. 
To limit the halo of the star so that it is not brighter than the planet, both turbulence 
and defects of the instrument have to be corrected with unprecedented precision, so a 
high-performance adaptive optics system is needed. The latest developments in adaptive 
optics (the first light from NAOS, a system installed by ESO on the VLT, was received 
in 2001) suggest that a system for direct detection of extrasolar planets may be available 
on 8-10 meter ground-based telescopes within 5 years. 

Several ground-based projects are currently being developed or under study. Most 
of them rely on interferometers to image planets. There are also several space based 
projects, among which are Darwin (studied by ESA) and Terrestrial Planet Finder 
(studied by NASA). These space interferometers would be able to detect Earth- 
like planets in habitable zones and to perform spectroscopic observations of their 
atmosphere. Given the ambitious nature of both projects, NASA and ESA may 
collaborate on the final mission. 
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Figure 1. Distribution of projected mass MpSini (in Mj, left panel) and semi-major 
axis a (in AU, logarithmic scale, right panel) for the known extrasolar planets. 



2.2. Statistical properties of observed extrasolar planets 

So far, about 5% of the stars monitored in radial velocity surveys have been found to 
harbor at least one planet. This is a lower limit of course, as low-mass and/or long- 
period objects cannot be detected with the current techniques. 

2.2.1. Mass and separation distributions: 

The histograms of the projected mass M^sini and semi-major axis a of the known 
extrasolar planets are shown in Figure 1. 

The hghtest object has a mass of 0.042 Mj or 13.3 Earth masses, which is very 
close to the mass of Uranus. Although the radial velocity surveys are biased towards 
massive planets, 37% of the planets detected so far have a mass smaller than 1 Mj and 
only 13% have a mass larger than 5 Mj. Within a few AU from the central star, there 
is a deficit of massive planets. 

The distribution of a shows that there is a significant number of planets (12%) with 
a between 0.035 and 0.05 AU, i.e. an orbital period between 2.5 and 4 days (the so-called 
"hot Jupiters"). Only 4 objects have a smaller period, which suggests that there is a 
pile-up at a rather well-defined separation (^^ 0.4 AU) or period 3 days). Only very 
few objects with a mass larger than that of Uranus seem to be able to "leak" through 
what looks like a barrier. Since observations are biased towards short-period planets, as 
noted by, e.g., Marcy et al. (2005), there is clear rise of the distribution with separation. 
Beyond ~ 0.4 AU, the data appears to be subject to y/N noise with dependence on 
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Figure 2. Projected mass Mp sin i (in Mj) vs. semi-major axis a (in AU) on 

logarithmic scales for the known extrasolar planets. Planets in multiple systems are 
indicated by open symbols. The dotted and dashed lines represent the detection limit 
for planets on circular orbits and an accuracy of 3 and 10 m s~^, respectively (eq. [1] 
with e = 0). 

binning, an approximate smooth representation for the number of planets N{a) with 
semi-major axis less than a in AU is dN / d(\oga) = 50 log a + 95. According to this, 
the number expected within an interval of 0.2 in log a is ~ 9 or 19 at logo = — 1 or 0, 
respectively. 

Figure 2 shows the projected mass MpSini vs. semi-major axis a. Planets in 
multiple systems are indicated by open symbols. There is a deficit of sub-Jupiter mass 
planets at separations larger than ~ 0.5 AU. More precisely, there is no planet with 
MpSini < 0.7 Mj at a > 0.5 AU. Since about 10% of the planets have a mass lower 
than 0.7 Mj, the absence of such planets at separations larger than 0.5 AU, where 57% 
of the planets are, is statistically significant. The dotted and dashed lines in figure 2 
indicate the detection limit for planets on circular orbits and an accuracy of 3 and 
10 m s^^, respectively. It is not clear that the deficit of low-mass planets at larger a 
is the result of observational biases (see also Udry et al. 2003). There is also a deficit 
of high mass planets at separations smaller than ~ 0.3 AU. More precisely, there is no 
planet with Mp sin i >4Mjata<0.29 AU. As observations favor massive planets 
with short-periods, this deficit is not due to an observational bias. It has further been 
pointed out that planets more massive than 2 Mj at a < 4 AU are actually members of 
binary star systems (Zucker & Mazeh 2002). 



2.2.2. Eccentricity distribution: 
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Figure 3. Eccentricity e i)s. semi-major axis a (in AU) on a logarithmic scale {left 
panel) and projected mass MpSin?' (in M.j) on a logarithmic scale {right panel) for the 
known extrasolar planets. Planets in multiple systems are indicated by open symbols. 



Figure 3 shows the eccentricity e vs. semi-major axis a and projected mass Mp sini. 
Planets with a < 0.05 AU have almost circular orbits, which is consistent with tidal 
circularization. At larger separations, the eccentricity varies from up to 0.93. This is in 
sharp contrast with the planets of our solar system, which are on nearly circular orbits. 
The upper limit of 0.93 may be special, as it may be due to the secular perturbation 
of the planet HD 80606 b by a stellar companion on an inclined orbit (Wu & Murray 
2003). Nonetheless, eccentricities in the range 0-0.8 are common. Note that higher 
eccentricities seem to be obtained for more distant planets. The eccentricity distribution 
of planets in systems of multiple planets is similar to that of planets in single systems. 

Figure 3 shows that higher eccentricities are obtained for more massive planets. 
This is consistent with the trend noted above of higher eccentricities at larger separations 
together with the fact that more massive planets are at larger separations. But note 
that the observed distributions could be produced either by physical processes tending to 
produce higher eccentricities for higher masses or higher eccentricities at larger orbital 
separations or a combination of the two. A suggested process of the former type is 
eccentricity generation through disk planet interaction ( Artymowicz 1992) and one 
of the latter type is orbital circularization close to the central star. Here again, the 
eccentricity distribution of planets in systems of multiple planets is similar to that of 
planets in single systems. 

Halbwachs et al (2005) suggest that the differences between the eccentricity 
distributions of planets and stellar binaries is an indication that these two classes of 
objects have different formation mechanisms. 
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Figure 4. Distribution of metallicity [Fe/H] for the stars harboring planets. 

2.2.3. Metallicity of stars harboring planets: 

Figure 4 shows the metalUcity distribution of the stars around which planets have 
been found. The metaUicity of most of the field dwarfs with no known companion lies 
in the range -0.3-0.2 (Santos et al 2001). The planet host stars are therefore metal 
rich, and the planetary frequency is clearly increasing with metallicity. This metallicity 
excess appears to have a "primordial" origin (Pinsonneault et al 2001; Santos et al 2001; 
Santos et al 2002; Sadakane et al 2002), rather than being due to the infall of planetary 
material onto the star which might result in observed heavy element polution for stars 
in the mass range 1-1.3 M0 (e.g. Sandquist et al. 1998, 2002). The favored explanation 
for the high metallicity of planet host stars is therefore that planets form more easily in 
a metal rich environment. 

Note that no short-period planets were found in the metal poor globular cluster 
47 Tucanae in which about 34,000 stars were monitored by HST. Simulations had shown 
that about 17 hot Jupiters should have been detected by photometric transit signals if 
the frequency of these objects were the same in the solar neighborhood and in 47 Tucanae 
(Gilliland et al 2000). 

Figure 5 shows the stellar metallicity vs. a, Mp sin i and e. Sozzetti (2004) has found 
that there is a weak correlation (at the 2 to 3a level) between the stellar metallicity and 
the orbital period of the system. This correlation becomes stronger when only single 
stars with one detected planet are considered. 



2.2.4. Planets in multiple stellar systems: 
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Figure 5. Stellar metallicity [Fe/H] vs. semi-major axis a (in AU) in logarithmic 
scale {left panel), projected mass Mp sin i (in Mj) in logarithmic scale {middle panel) 
and eccentricity e {right panel) for the known extrasolar planets. Planets in multiple 
systems are indicated by open symbols. 



As noted above, the most massive short-period planets are all members of stellar 
binaries (Zucker and Mazeh 2002). It also appears that planets orbiting a component of 
a multiple star system tend to have a very low eccentricity when their period is shorter 
than about 40 days, i.e. when a is smaller than about 0.2 AU (Eggenberger et al 2004). 

2.2.5. Free floating planets: 

Recently, a population of about 30 very young (a few Myr) free-floating objects with 
planetary masses (from a few to 13 Mj) has been found in the Sigma Orionis open clus- 
ter (Zapatero Osorio et al 2000) and in the Trapezium (Lucas and Roche 2000). The 
lightest of these objects is S Ori 70 (Zapatero Osorio et al 2002). Its mass is derived 
from evolutionary models and is estimated to be 3 Mj. However, Burgasser et al (2004) 
claim that this object may actually not be a member of the a Orionis cluster, in which 
case its mass may be much higher. 

Bejar et al (2001) have pointed out that if the distribution of stellar and substellar 
masses in a Orionis is representative of the Galactic disk, older and much lower 
luminosity free-floating objects with masses down to a few Mj should be abundant 
in the solar vicinity, with a density similar to M-type stars. 

The discovery of these free-floating objects has brought some confusion in the 
definition of planets, as they are not orbiting a star and they probably have not formed 
through the standard core-accretion scenario (see below) . 

3. Theories of planet formation 

Terrestrial planets are believed to be formed via solid body accretion of km-sized objects, 
which themselves are produced as a result of the sedimentation and collisional growth of 
dust grains in the protoplanetary disk (see Lissauer 1993 and references therein). This 
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theory comes from the idea already proposed in the 19th century that the Earth and 
the other terrestrial planets were formed from meteoritic material. At the beginning of 
the 20th century, this idea was supported by the geologist T. C. Chamberlin who called 
"planetesimals" the solid objects that had accumulated to form the planets. But it is 
only in the 1960s that the theory was developed in a quantitative way by V. Safronov 
(1969), who calculated in details the different stages of terrestrial planet formation. 
This pioneering work, originally written in Russian, was published in English in 1972 
under the title "Evolution of the Protoplanetary Cloud" . Since then, the theory has 
been further developed, in particular by G. W. Wetherill and collaborators who have 
used intensive numerical simulations to study planetesimal accumulation. 

As far as giant planets are concerned, two theories have been proposed. According 
to the first theory, giant planets are formed through collapse and fragmentation of 
protostellar disks. This scenario is related to that proposed by Laplace in 1796. 
In its modern version, it has been developed by G. Kuiper in the 1950s and by 
A. G. W. Cameron in the 1960s and 1970s. It is still being studied today. The second 
theory, proposed by Cameron in 1973, is called the "core accretion model" . A sohd core 
is first assembled in the same way as terrestrial planets. Once the core becomes massive 
enough to gravitationally bind the gas in which it is embedded (typically at around a 
tenth of an Earth mass (see section 3.2.2 below), a gaseous envelope begins to form 
around the core. 

3.1. Terrestrial planet and planetary core formation 

In the formation of terrestrial planets or planetary cores, several stages can be 
considered. Firstablc, micron sized dust grains that are present when the disk forms 
sediment towards the disk midplane. These grains are subject to both the drag force 
exerted by the gas and the gravitational force from the central star. Accumulation of 
grains during this first stage leads to cm/m sized particles. It is still not clear how these 
objects further grow to become 0.1-1 km sized planetesimals. It is likely that growth 
occurs through collisions and accumulation. The planetesimals are massive enough 
that mutual gravitational interactions have to be taken into account in computing their 
evolution. Further collisions and accumulation lead to terrestrial planets or the core of 
giant planets. These different processes are described in details below. 

3.1.1. Growth from dust to cm/m sized particles: 

When the disk forms out of the collapse of a molecular cloud, it contains dust grains 
which size is that of the interstellar grains, i.e. 0.1-1 /im. Under the action of the 
gravitational force exerted by the central star and the drag force exerted by the gas, 
when there is no turbulence present, these grains tend to undergo damped oscillations 
around the disk midplane, toward which they eventually sediment. Simultaneously, they 
drift toward the disc center. This is because the coupling with the gas forces them to 
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rotate at a slightly sub-Keplerian velocity. Since they are not supported by pressure 
forces, they consequently cannot stay on a fixed orbit. 

For typical disk parameters, the gas mean free path around the disc midplane is 
~ 1 m. This is very large compared to the size of the grains we are interested in here. 
In this regime, the drag force exerted by the gas on the grains is ~ —■nr'^pcs (v — v^) 
(Epstein law; e.g., Weidenschilling 1977), where is the radius of the grains, which are 
assumed to be spherical, Cg is the thermal velocity of the gas (i.e. the sound speed), 
p is the mass density of gas in the disk, is the velocity of the grains and v that of 
the gas. This formula is vahd provided |v — v^l ^ c^, which is the case here as can be 
checked a posteriori. We define the time: 



Trrjpcs '' 



with rud being the mass of a dust grain. This is the timescale it takes for the velocity 
of a grain subject only to the drag force to decrease by a factor e. For ra — A*m, 
p ~ 10~^° g cm~^, Cg ~ 10^ cm s~^ (which corresponds to a temperature of about 
700 K), which are typical conditions at about 1 AU from the central star, and assuming 
the mass density of the grains to be 3 g cm^'^ (which is comparable to the mass density 
of the Earth's crust and that of the asteroids), we get Tg ~ 10 s. 

For the gas, we write v = {vriV^^iVz) in cylindrical coordinates and, for the dust, 
Vd = {vr^d) v^^d, ^z,d) ■ We assume the motion of the grains to be almost Keplerian with a 
small drift velocity, i.e. v^^d ~ vk and \vr,d\ ^ k^j.dL where vk is the Keplerian velocity. 
In this situation in a steady state, the dust to a first approximation comoves with the 
gas so that v^^d — v^- This means that accelerations due to radial pressure gradients 
acting on the gas have to be balanced by the accelerations due to drag acting on the 
dust which occur because of the radial drift velocity which will be given by: 

Te r ^ p or 

As expected, if the pressure increases inwards, Vr^d < 0. Also Vr^d ^ as — > 0, which 
means that very small grains are "carried" by the gas. With \dP/dr\ ~ P/r, cl P/p 
and H/r ~ Cs/vk, where H is the disk semi-thickness, we get Vr,d ~ —Cs{H/r)Te^K-i 
where Vt^ is the Keplerian angular velocity. With the values of the parameters used 
above, this gives v,.,d on the order of a few 10^^ cm at r = 1 AU. 

If the grains were not subject to the gas drag force, they would oscillate around the 
disk midplane with the frequency VLk- Since 2t{ /VLk ^ Te, the oscillations are suppressed 
and the grains sediment towards the disk midplane. As for a highly damped harmonic 
oscillator, the characteristic timescale for sedimentation is 

1 3pcs 



nj^Te ^^KpdTd 

Here pd is the mass density of the grains. For the values of the parameters used above, 
Tg — 10^ years at 1 AU. This is small compared to the planet formation timescale, as 
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will be seen below. We have supposed that md and are constant, so that given 
above is actually a maximum. Indeed, if the grains grow through collisions as they fall 
toward the midplane, the drag force due to the gas decreases and the sedimentation is 
accelerated. Note that during Tg, the grains drift radially over a distance |fr,d|Ts, which is 
on the order of 10~^ AU at 1 AU. This is very small compared to the disk semithickness 
(~ 0.1 AU), which justifies the neglect of Vr,d in the discussion of the vertical motion. 
However, grains grow through collisions during their descent. Let us suppose that they 
"absorb" all other grains on their way to the midplane. When a grain moves a distance 
\dz\, its mass then increases by dma — Trr'^\dz\ps, where Ps if the mass density of solids 
(dust) in the disk. This corresponds to an increase of its radius of drd = ps\dz\/ {Apd). 
Therefore, for the values of the parameters used above, a grain initially at the altitude 
H would have a size ~ mm when it settles toward the midplane at a few AU. Its radial 
drift velocity then cannot be neglected. 

Numerical simulations in laminar disks (Weidenschilling 1980, Nakagawa et al 1981) 
show that grains in the process of settling toward the midplane undergo a significant 
radial drift which promotes collisions. This enables them to grow further up to 0.1-1 m, 
and sedimentation is then found to occur in a few 10^ years (sedimentation being faster 
for larger grains, as Tg oc l/r^). Note however that these simulations assume that the 
grains stick perfectly to each other, so that they agglomerate when they collide. 

Even though the growth of the grains during their sedimentation is not yet fully 
modeled by numerical simulations, it is generally accepted that this process leads to 
cm/m sized particles around the (laminar) disk midplane in less than 10*^ years. 

Calculations including turbulence have been performed by Weidenschilling (1984). 
He found that the coagulation rate is initially much greater than in a non turbulent 
disk. Aggregates quickly reach sizes ~ 0.1-1 cm, but erosion and breakup in collisions 
prevent growth of larger bodies. These aggregates are too small to settle to the plane 
of the disk in the presence of turbulence. 

3.1.2. Growth from cm/m sized particles to 0.1/1 km sized planetesimals: 

Before examining the processes that could lead to further growth of the solid par- 
ticles, let us discuss the timescale on which such processes may have to operate. 

Maximum radial drift: When the gas mean free path I is no longer large compared 
to the size of the grains r^, the Epstein law cannot be used to calculate the drag force 
between the gas and the grains. In general, the drag force acting on a single grain is 
given by — / (v^ — v), where / is a coefficient that depends on the characteristics of the 
grains and the gas. In steady state, the horizontal equations of motion for the grains 
can therefore be written: 




(6) 



r 
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Vr,d = - —^V<p, (7) 

ma 

where Av^ — v^^d — v^- Note that: 

In equations (7) we have used ^ \''^tp,d\ l^'^^'l ^ '^v'- We can further ehminate 
f/rrid to get a quadratic equation for Av^. It is straightforward to show that this 
equation has a solution only if 
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^ 2d_ f]_dP\ ^ ^ 
pr dr r dr \pr dr I 2 



vL > 0- (9) 



The first two terms in bracket are negligible compared to the last one since Cg/ {rflx) ^ 
H/r <^ 1. This condition can therefore be written \vr,d\ < \v:f — vk\ ~ c^/(^^x), where 
again we have used v^/vk — i- At 1 AU, this gives a maximum value of \vr,d\ on the 
order of 10^ cm s~^, independent of r^. Prom subsection 3.1.1, \vr,d\ ~ 10~^ cm s~^ 
when rd ~ l//m <S i.e. when the Epstein law for drag still holds, much smaller than 
the maximum speed derived here. As \vr,d\ r^, the maximum speed is reached when 
rd is at least as large as / ~ 1 m. The inflow time corresponding to the maximum speed 
is very short, being less than 100 years at 1 AU. 

Radial drift when rd^ I: In the regime 3> the drag force acting on a single 
grain is ~ —CoT^r'^p |vd — v| (v^ — v) , where Co is the drag coefficient, being a number 
of order unity (Stokes law; e.g., Weidenschilling 1977). In steady state, the radial 
equation of motion for the solid particles leads to v^p^d — if we assume \vr^d\ ^ \''^<p,d\- 
This means that the angular velocity of large particles around the star is not affected 
by the gas. The azimuthal equation of motion (7), after using equation (8), then gives: 



'Jr,d 



Corl (dP\ 



— 



2mdP^K \ ^'^ / 



dP 



dr 



(10) 



Important results follow from the dependence of the migration speed on the pressure 
gradient. In particular, if the pressure increases inwards, inward drift occurs. Similarly, 
an outwardly increasing pressure results in outward drift. Thus, particles will 
accumulate at pressure maxima and be depleted at pressure minima. We comment 
that this is a general feature of a non magnetic gas and dust mixture that arises because 
the forces acting on the two components arc the same at pressure extrema so that they 
can comovc there without drag. To estimate the typical drift speed magnitude in the 
general case, we set dP/dr = nP/r. Then equation (10) gives: 

_ 3Cgn|n|S eg 
'^^'''""^ 16 r^V ^ ^ 

where S is the surface density of the gas disk. 

For Ce> = 0.2, \n\ = 3, pd = 1 g cm~^, = 10 m and S = 10^ g cm~^, wc get 
\'^r,d\ = 0.1cs[cs/(rr2x)]^. The radial flow time scale is thus {5/n){rflK/cs)^Porb, with 
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Porb being the local orbital period. For Cs/{rVtK) = 0.07, this is characteristically 
~ 5000 orbits, being more than an order of magnitude longer than the timescale 
corresponding to the maximum drift speed. Therefore, the maximum speed is reached 
when ~ / ~ 1 m. 

Although there is a dependence on the local state variables, the drift timescale for meter 
sized bodies at 1 AU is potentially as short as ~ 100 years. Therefore, either the process 
by which sub-meter sized bodies grow to become planetesimals has to be very fast, or 
possibly special flow features, like pressure extrema or vortices, play an important role 
(see below). 

Note however that the calculation presented above is for a particle which is isolated. 
In a laminar disk in which a dust layer forms around the midplane, the drift velocity 
cannot be calculated by assuming that the gas velocity is given by equation (8). In the 
dust layer, the dust density is probably large enough that the gas is carried around at 
a velocity closer to the Keplerian value. The relative velocity of the solid particles with 
respect to the gas is therefore smaller than what we have assumed. It is actually the 
collective drag on the surface of the layer which is responsible for orbital decay. The 
hfetime of the dust layer against orbital decay at 1 AU is estimated to be ~ 10^ years 
(Goldreich and Ward 1973) so that the problem still remains. 

The process by which 0.1/1 km sized bodies form in the dust layer in a laminar disk 
is still unclear. Safronov (1969) and Goldreich and Ward (1973) independently pointed 
out that the dust may concentrate sufficiently in the layer to undergo gravitational 
instability. This occurs when the Toomrc parameter Q = Qcd/iirGHd) (with Cd and 
T,d being the velocity dispersion and surface mass density, respectively, of the dust) 
approaches unity. In a typical disk, this requires the mass density of dust to be 
~ 10~^ g cm~^. Fragmentation of this layer then produces planetesimals which size 
is on the order of 1 km. The advantage of this scenario is that it operates very fast, on 
a timescale on the order of the orbital timescale. 

However, it was pointed out by Weidenschilling (1980) that, if the dust density 
reaches such high values, then it dominates the disk in the midplane, i.e. it carries the 
gas at a Keplerian velocity around the star. As mentionncd above, because there is 
no pressure associated with the dust, it rotates at a different rate than the gas around 
the layer. In these conditions, the velocity shear between the dust layer and the gas 
around it gives rise to Kelvin-Helmholtz instabilities that result in turbulence, even 
in a disk which is laminar to begin with. This leads to mixing of the dust so that 
the density in the layer is decreased and gravitational instabilities are prevented. This 
result was confirmed by detailed calculations by Cuzzi et al (1993) who showed that 
rapid accretion of planetesimals by gravitationally unstable fragmentation on an orbital 
timescale is unlikely to occur until objects have already accreted by some other process 
to the mass of the largest known meteorite samples, if at all. Note that if the disk is 
turbulent to begin with, settling is inhibited for particles below a certain size ~ 1 cm, 
as already mentioned above. 

The instability scenario has recently been revisited. Sekiya (1998) showed that 
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gravitational instabilities can occur if the dust to gas ratio is much larger than the cosmic 
value. Youdin and Shu (2002) have discussed mechanisms for increasing this ratio. In 
particular, they showed that gas drag in a laminar disk causes global redistribution and 
concentration of small solids as they drift inwards. This is because the drift velocity 
of these particles decreases as they move in to smaller radii. Youdin and Shu (2002) 
conclude that the dust to gas ratio can become high enough for gravitational instabihties 
to develop in less than a few miUion years. However, WeidenschiUing (2003) pointed 
out their arguements are based on motions obtained for independent particles and 
claimed that collective motion due to turbulent stress on the particle layer acts to 
inhibit concentration of particles and may prevent gravitational instability. 

Recently, Goodman and Pindor (2000) discovered a secular instability that may 
lead to rapid formation of planetesimals in the dust layer. This instability is not driven 
by self-gravity but by drag. If the drag exerted by the (turbulent) gas on the surface 
layer is collective, then the drift speed of the dust layer varies inversely with its surface 
density. A slight increase of the mass of some annulus would therefore result in a 
decrease of its drift velocity, leading to the accumulation of the dust inflowing from 
adjacent larger radii. The instability, which grows exponentially in a frame moving at 
the unperturbed drift speed, leads to the formation of overdense rings on an orbital 
timescale. These rings ftuthcr fragment to form planetesimals of size ~ 10 km at 1 AU. 
However, turbulent stresses acting on the particle layer may act to inhibit concentration 
of particles in this case also. 

If none of the mechanisms described above operate, one has to rely on sticky particle 
collisions (coagulation) to obtain planetesimals. Cuzzi et al (1993) argue that 10-100 km 
sized objects can be formed through this process in the layer in about a million years (also 
see below). Note that coagulation and growth can be made more efficient if particles 
can be trapped in vortices or overdense regions of the disk. Planetesimal formation in 
vortices has been studied by, e.g.. Barge and Sommeria (1995), Tanga et al (1996) and 
Klahr and Henning (1997) who found that growth was indeed much faster. But note 
that Cuzzi et al. (2001) make the point that dust density enhancement may only be 
significant for a narrow particle size range. Furthermore, it is not clear that vortices 
can live long enough in disks for any dust densification process to be efficient. More 
recently, trapping of particles in the spiral arms of a gravitationally unstable disk has 
been examined. Numerical simulations of 1-10 m sized objects in marginally stable, 
self-gravitating disks by Rice et al (2004) show that the drag force causes the solid 
objects to drift towards the peaks of the spiral arms where the density and pressure are 
highest, in accord with the behavior predicted by using equation (10) (self-gravity was 
not taken into account in this equation, but it would not be expected to modify the 
situation significantly as it acts equally on dust and gas). It is further speculated that 
the density enhancements may be sufficient for the growth of planetesimals through 
direct gravitational collapse. Similarly, Durisen et al (2005) suggest that the dense 
rings that appear in their simulations of gravitationally unstable disks are conducive to 
accelerated growth of planetesimals. 
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However, there remains an issue as to the hfetime of these features in an unstable 
self-gravitating disk. In future years, more may be learnt about the feasibility 
of the processes discussed here through more extensive three dimensional numerical 
simulations. 

3.1.3. Accretion of a planetary core from planetesimals: 

Once planetesimals of radius ~ 10 km have formed, they proceed to evolve by pro- 
cesses of velocity dispersion growth due to mutual gravitational interaction, gas drag, 
disk planet interaction and coUisional accretion of mass which results in the accumula- 
tion of massive solid cores. 

It is believed that the process starts with some planetesimals entering a regime of 
runaway growth and forming a population of embryo cores. Once these are massive 
enough to significantly increase the velocity dispersion equilibrium of the smaller 
planetesimals, runaway growth ceases and the system enters an oligarchic growth phase 
in which cores separated by 15 or so Hill radii accrete at comparable rates. Throughout 
these phases, gravitational focusing plays an important part in making the accretion 
rates rapid enough that cores can enter the Earth mass range within the lifetime of the 
disk on a scale of 5 AU. 

After this, the cores can become isolated and processes like gap formation may 
come into play to inhibit further accretion. The mobility of the cores as a result of disk 
protoplanet interaction in the presence of turbulence may be important. Note that, if 
the disk is laminar, rapid inward type I migration due to disk protoplanet interaction 
poses a serious threat to core survival. Thus the final accumulation of cores up to the 5- 
15 range is still poorly understood. Nonetheless, if favourable circumstances exist 
(see discusion and references in sections 4.2.1 and 4.2.2 below), this may be possible 
within the lifetime of the gas disk and thus giant planets can form. 

A population of cores may be left behind to undergo gas free accumulation when 
the gas disk has disappeared. The formation time of the Earth through this process 
is expected to be ~ 2 x 10^ yr. In the outer solar system, gravitational scattering of 
the cores can produce escapers leaving a residue to accumulate the remaining small 
planetesimals and form Uranus and Neptune in a few billion years. 

We now go on to discuss these processes in more detail. In doing so, we follow 
many authors and adopt a so-called standard, or minimum mass, solar nebula model. 
This has a gas surface density profile E oc r~^/^ and 2 Jupiter masses within 5.2 AU 
(Hayashi 1981). This is an estimate of the spread out minimum mass required to form 
the planets in the solar system. We further adopt a solid or condensate surface density 
of 1% of the gas surface density. The characteristic lifetime of the gas disk is between 
one and ten million years and giant planets have to be formed on this timescale in the 
solid core followed by gas accretion scenario. 

The gas disk is also believed to undergo outward angular momentum transport 
allowing accretion onto the central star in the manner of an accretion disk. Associated 
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with this is an effective kinematic viscosity, u, the origin of which is most likely associa- 
ted with turbulence driven by the magnetorotational instability (Balbus and Hawley 
1991, 1998). In the standard model, we adopt z//(r^fi) = a{H/rYVl = 10"^, with a 
being the well known Shakura & Sunyaev (1973) dimensionless viscosity parameter. 

• Gravitational focusing and runaway accretion: 

We first show that direct accumulation of planetesimals cannot lead to the formation of 
Earth sized objects within the disk lifetime. Consider the simplest possible accumulation 
scenario where a core of mass m„ grows through the accretion of sticky solid particles 
of mass rUp which have number density Up, and thus mass density Pp = UpTUp, through 
direct impacts. We suppose that the accreting particles are decoupled from any gas 
present, so that their local root mean square velocity dispersion, Vp, is related to their 
estimated vertical semithickness. Hp, through Hp = Vp/fl and their surface density is 
Yjp = IppHp. 

Neglecting any effects due to gravitational focusing of the colliding masses, increases 
at a rate given by: 

= TTO PpVp, (12) 

where a is the core radius. Thus the mass accumulation timescale is given by: 

1 ^ dmg 1 ^ 3SpQ ^^^^ 
tacc 8apa ' 

where the density of the accumulating core is Pa- This gives: 

Adopting Pa — i gm cm~^ and Sp/E = 0.01, we find: 

tacc = 40 ( ^ ^ (^f\r. (15) 

\1 g cm-2 J 1 cm VI AUy ^ ^ ' 

This indicates that, for disk models like the minimum mass solar nebula model with 
E ~ 200 g cm~^ at 5 AU, straightforward accumulation processes operating with ef- 
ficient sticking could conceivably form objects of up to 10-100 km in size at 5 AU, 
within the lifetime of protoplanetary disks of a few million years. However, to produce 
objects of significantly larger size, an increase of the accretion cross section resulting 
from gravitational focusing, which can occur at low relative impact velocities, needs to 
be invoked. 

To take into account gravitational focusing, we have to include the gravitational 
attraction of the core on the colliding particles. This increases the impact parameter 



for a collision from the radius a to a^s — a 2Gma/ {av, 



1/2 



The effect of this 
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gravitational focusing is to increase the effective cross section for impacting collisions 
accordingly. Thus equation (13) is modified to become: 



The enhancement factor given by the term in brackets on the right hand side of 
equation (16) is large when Vp is very much smaller than the escape velocity from 
the core. For one Earth mass with a mean density of unity at 5 AU immersed in a 
planetesimal swarm with Vp = 0.01rr2, the enhancement factor is close to 4250. 
Thus, under these conditions, forming cores in the Earth mass range in disks with 
properties similar to those of the minimum mass solar nebula becomes a possibility. 
Whether the relative impact speeds are at an appropriate level depends on the velocity 
dispersion equilibrium attained by the distribution of cores and impacting particles or 
planetesimals. This we shall consider below. 

At this point, we comment that when gravitational focusing dominates, a phe- 
nomenon known as runaway accretion is possible (e.g., Safronov 1969; Wetherill & 
Stewart 1989). This occurs when the accretion rate increases rapidly enough with mass 
such that, if a particular core has a mass slightly exceeding that of other cores, it will 
grow increasingly faster and so run away from its neighbors in mass. 
For example, consider the situation when Vp is fixed. Then, given a oc m]/^, from 
equation (16) we obtain dm^/dt ex m^, with k — 4/3. When, as in this case, k > 1, 
rria — > oo in a finite time, such that most of the evolution time is spent with rria in the 
neighborhood of the smallest mass followed by rapid growth to large values. Hence the 
term 'runaway'. In this situation, ultimately a single mass dominates. More generally, 
runaway accretion requires Vp oc m^' with k' < 1/6. However, note that accumulation 
with strong gravitational focusing may occur without runaway. Then most time is not 
necessarily spent with rria close to the smallest mass. 

Determination of the root mean square velocity dispersion Vp for the impacting 
planetesimals requires consideration of the balance between growth due to gravitational 
scattering and damping due to interaction with the gas which we now consider. 

• Analysis of Gravitational Scattering: 

When the planetesimals which form the building blocks from which the core must accrete 
have a single characteristic mass rria ~ 10^^ g, we require there to be ~ 10^° such objects 
that can be accumulated to form a body in the Earth mass range. Such a large number 
is best treated by statistical means. A system of many bodies interacting under their 
mutual gravitation can be described by the Fokker-Planck equation which we write in 
the form: 




(16) 




(17) 
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Here /„ denotes the phase space number density of bodies with mass rria- The operator 
giving rise to evolution due to mutual interactions among the planetesimals which lead 
to gravitational scattering is Fcou. The operator causing evolution due to interaction 
with the gaseous disk is Fgas- The latter combines the effects of gas drag, eccentricity 
and inclination damping as well as orbital migration. The derivative operator is taken 
following a particle orbit under the gravitational force due to the central mass. Thus: 

D d d „^ d 

where the central potential at a distance r from the central mass is $ = —GM^/r. 
The particle position and velocity vectors measured with respect to an arbitrary 
coordinate system are denoted by r and v, respectively. 

We consider the situation when the planetesimals have a dispersion velocity that is very 
small compared to the orbital velocity. In this case, the position of a particular plane- 

tesimal does not deviate much from a particular radius on orbital timescales, implying 
that a local description of the interactions is possible. Such a description can be given 
in the context of a local shearing box (Goldreich & Lynden-Bell 1965). 

— Local shearing sheet approximation: Distribution function for the planetesimal 
swarm: In the shearing box, we consider a uniformly rotating local Cartesian coordinate 
system for which the origin, located at some point of interest in circular orbit with 
radius/semimajor axis r, corotates with the Keplerian angular velocity Q. The x-axis 
points radially outwards, the y-axis points in the azimuthal direction in the direction 
of rotation while the ^-axis points in the vertical direction. A linear expansion for the 
combined acceleration due to gravity and the centrifugal force is used such that: 

_ V - ^qV^) = (3n^x, 0, -n^z) . (19) 

Then for an axisymmetric disk with no dependence on y: 

D/a dfa , dfa , df„ 



Dt dt dx dz 

_2nv,^ + (3nx + 2vy)n^-n^z^. (20) 

OVy OVx OVz 

The Liouville equation D/a/Dt = does not have a solution corresponding to an 
isotropic Gaussian but does have one for which takes the form of an anisotropic 
Gaussian given by: 

Here the root mean square velocity dispersions (cr^:, cTj,, cr^) are constant and such that 
= (T2:/2. There is no constraint on az- The velocity Uy = Vy + 'iVtx/2 is measured 
relative to the local circular velocity, Vy = —3flx/2, and we shall adopt local relative 
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velocity vectors v = {v^, Uy, v^). The constant Cq, is related to the spatial number density 
in the midplane, ria, through: 



expected to accrete from a planetesimal swarm with a distribution of masses that itself 
evolves under gravitational scattering and accretion. However, following Wetherill & 
Stewart (1989) and Ida & Makino (1993), we simplify matters by assuming that there 
are only two populations present. One corresponds to the growing cores and the other 
to the planetesimal swarm from which they accrete. The latter distribution consists 
of a large number of objects with a small fixed mass, rua- The former consists of a 
much smaller number of cores of larger mass, m^, which increases as the system evolves. 
We adopt rria = 10^^ g. This is the mass that can be accumulated within the disk 
lifetime without the need for gravitational focusing as indicated above. This is also the 
mass that emerges from the assumption that planctcsimals are able to form through the 
gravitational instability in a dust layer (Goldreich and Ward 1973). 
The evolution of each component may be described by a Gaussian such as that given 
by equation (21) with velocity dispersions and number density that slowly evolve under 
the action of collisions. 

— System evolution through gravitational scattering: Encounters between planetesimals 
that occur without direct physical impacts tend to convert kinetic energy from shear 
into that of random motions, in the same way as the action of viscosity converts energy 
from shear into thermal motions in a gaseous disk. For the planetesimal swarm, effects 
due to evolution resulting from gravitational scattering occur on a timescale that is 
much longer than orbital. 

To model this, we use the form of Fcou given by Binney & Tremaine (1987). This neglects 
rotation about the central mass, which should be a reasonable approximation as long as 
the timescale associated with a binary encounter is short compared to This in turn 
requires that cr^/Q > rn, where rn — [?7^-y/(3M^)]^/^r is the Hill radius appropriate to 

the largest characteristic mass, m^, involved in the encounter. 

Following Binney & Tremaine (1987), we write using the summation convention for 
repeated indices: 



We consider two kinds of interaction. The first, which is important during the early 
stages of the core mass build up, arises when objects in the distribution with small 
mass, rua, interact with each other. Then 




In general, a growing core is 



(22) 




(23) 



Ai = 47rG^ ln(A„)m; 



2 9h 



(24) 
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and 

A, = 4.GMn(A„)mJ^, (25) 

with 

g = j /,(v')|v - v'l dV and h = I dV. (26) 

We take = Sa^Ha/ {AGma) , giving the ratio of maximum to minimum impact 
parameters as disc semithickness to impact parameter for a typical deflection (e.g., 
Binney & Tremaine 1987). 

The second type of interaction we consider, which is important during the later stages 
of core mass build up, is between an object from the first distribution with small mass, 
rUa, and an object from the second distribution with large mass, m0, the number density 
of these being Because of their high inertia, these can be assumed to be in circular 
orbit and have zero velocity dispersion. In this limit: 

Ai = 0, (27) 

and 

Dij = AnG^rif:) ln(A;3)mJ (dijV^ - ViVj) /v\ (28) 
where v — |v|. Here = 3alHa/{4:Gmp). 



— Growth of the velocity dispersion: The effect of gravitational scattering is to cause the 
velocity dispersion to increase. This is most easily seen by formulating the Boltzmann 
Ti theorem for the problems we are considering. This states that for a single mass 
species interacting either with itself or another in fixed circular orbits with no velocity 
dispersion, 

n=- lfaHfa)d\ (29) 

increases monotonically with time. H, which can be related to the entropy, can remain 
constant only for an isotropic Gaussian which cannot be attained here because of the 
form of particle orbits in the central potential (see eq. [21]). Using the distribution 
function given by equation (21), and assuming that ria and the ratio of velocity 
dispersion components are independent of time (see Papaloizou & Larwood 2000), we 
obtain at the midplane (z = 0): 

dt ~ a^, dt ' ^ ^ 

Thus dH/dt > implies that the velocity dispersion must increase with time. 

We first evaluate dTi/dt for a single mass species interacting with itself. Using the 
expression (29) of 7i and rcoii(/a) given by equations (23), (24) and (25), we obtain: 
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+ 2.Gm„ln(A,)y^^-^^^dvdv. (31) 

When (Ty = (T2, the integral may be evaluated analytically for an anisotropic Gaussian 
specified by equation (21). The relaxation time defined through tu — (Tx/ {dcTx/dt) can 
then be calculated from equation (30) in the disk midplane: 



1 ^ 80FG'^m^ln(Ae,)ne 




2 + 
2- 



(32) 



Here we have used the relation ay — (7^/2 which is valid in the Keplerian case. 

In terms of the semithickness of the planetesimal disk Ha — cr^/^ and the planetesimal 

surface density Eq, = y/2nHamana, we have 



A„=,^)(^|, (33) 



' omHK)^(^X^, (34) 



and 

tR "—-'-"^ M2 \Ha 

where Mo — TrE^r^ gives an estimate of the mass within the protoplanetary disk 
comprising the planetesimal swarm that is contained within r. 

A number of effects may act to oppose the growth of the velocity dispersion induced 

by gravitational scattering and allow a quasi steady state distribution to be reached. 

These include (see below for additional details) gas drag, tidal interactions with the disk 

and physical collisions between planetesimals which, if very inelastic, will prevent the 

dispersion velocity from significantly exceeding the planetesimal escape velocity. 

To estimate likely relaxation times, we adopt M^, — 1 M©, rua — 10^^ g together with 

Ha/r = 3 X 10~^, which corresponds to a dispersion velocity comparable to the escape 

velocity, and = 10 M®. This gives ~ 5 x 10^ yr at r = 1 AU. 

But note that the relaxation time is very sensitive to the velocity dispersion, making 

the processes that determine it significant for in turn determining the timescale for core 

accumulation. 



We now evaluate dH/dt in the case when the gravitational interactions of the 
planetesimal swarm, labelled by a, with the system of cores with large masses in 
circular orbits with zero velocity dispersion, labelled by /3, become more important than 
interactions between its own members. Using equations (23), (27) and (28), conjointly 
with the expression (29) of H, we find: 

— = 2nG m,r^, ln(A,) j ^- d v. (35) 
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In this case, with identical assumptions about the ratios of the velocity dispersion 
components as in the self-interaction case, the integral is readily performed to give 
the relaxation time in the midplane through: 



In (2 + Vs) 



(36) 



tR,mid V^C^ 

Here = {3M^,H^)/{mpr^). When, as here, mp » m^, the disk semithickness 
associated with the heavy distribution, Hj3 — E/j/(-\/27rm/jn;3), is much smaller than 
that associated with the light distribution. Accordingly, we take an appropriate average 
of the relaxation rate by multiplying it by the ratio SHp/ {^/2Ha)j this being a measure 
of the proportion of time spent by a member of the light distribution interacting with 
the heavy distribution. Thus we obtain the final mean relaxation rate: 

^ln(2 + V3)-l 



1 _ 2AG^mpln{Ap)J:p 



tR V^H^al 



2 



(37) 



— Comparison of scattering by small and large masses: As an ensemble of planetesimals 
evolves, the direction of evolution is expected to be such that, because of runaway 
accretion, the distribution becomes dominated by an increasing number of larger masses. 
Prior to this, when the smaller masses dominate the dynamical relaxation, because their 
mass does not change, and their velocity dispersion equilibrium is determined by balance 
between dynamical relaxation and damping through interaction with the gas disk, the 
relaxation rate is independent of the magnitude of the larger mass mp. Accordingly, 
runaway accretion by these may occur (see above). 

The condition for the larger masses to dominate the relaxation is determined by equating 
the relaxation rates given by equations (32) and (37) and is, neglecting any variation in 
ln(A): 

3m/3S/5 > m^SQ. (38) 

Simulations (see, e.g., Ida & Makino 1993, Thommcs et al 2002) indicate that, once the 
larger masses start to dominate the relaxation, the system enters the orderly mode of 
oligarchic accretion in which the masses mjs remain in near circular orbits separated by 
~ 10/ times the Hill radius rn = [m^/(3M^)]^/^r, with / being a number of order 
unity, and grow by capturing the low mass material in their neighborhood. Then 
Tip — mp/ {20 fnrrH)- Using this, the condition (38) implies that for ohgarchic accretion 
we require: 

> 2.2 X {^f (j^^^^f M,. (39) 

This strongly suggests that the transition to oligarchic accretion occurs at quite an early 
stage in the core accumulation process. 



— Damping processes: 
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Gas Drag: We consider a planetesimal of mass rria, density pa and radius a moving 
relative to the disk gas with speed v. The equation of motion for v is: 

ma— = -na^Copv^, (40) 

where Cd, being a dimensionless number normally in the range 0.1-1, is the drag 

coefficient. This leads to a relative velocity damping time, tgd = v/{dv/dt), equivalent 
to: 

1 _{367iy/^ECDer ^^^^ 



where we have set v — erfi ~ with e corresponding to an orbital eccentricity for a 
planetesimal. In cgs units, this gives: 

Note that this time can be small for small masses, which is related to the need to evolve 
fast through the size range where decouphng from the gas first occurs to objects of 
~ 10^^ g (see section 3.1.2). 

Relative importance of gas drag and disk tides: Because of the dependence on 
protoplanet mass, gas drag is more effective for smaller masses while disk tides take 
over for larger masses. The cross over mass, which is also the mass for which the 
damping time is a maximum, can be estimated by equating the gas drag time given 
by equation (42) with the eccentricity decay time resulting from disk tides given by 
equation (78) (see below). This gives: 



Me ■ VO-05 J \ e J V 5.2AU 

Although there is significant sensitivity to protoplanet parameters, disk location, eccen- 
tricity and aspect ratio (going from 5.2 AU to 1 AU, other things being equal, increases 
it by an order of magnitude) the crossover mass is characteristically on the in the range 
lO^'^ — 10^^ earth masses. Thus gas drag will be the dominant eccentricity damping 
process when the planetesimals first form. However, if many larger mass cores form and 
then interact during the later stages of core accumulation, disk tides will dominate. 



• Velocity dispersion equilibrium: 



From the above estimates we anticipate the following evolutionary sequence: In the 
early phases, when the core masses are comparable to the planetesimal mass, the ba- 
lance is between relaxation by self-interaction among the small masses against gas drag. 
In this case the velocity dispersion depends only on the small unchanging mass nia- A 
few masses can then undergo runaway accretion growing to large values mp. When this 
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is high enough, the larger masses dominate the gravitational scattering of the planetes- 
imals which is again balanced against gas drag. The velocity dispersion in this phase 
increases with m^j and the core accumulation moves into the more orderly oligarchic 
regime (Ida & Makino 1993). 



— Scattering due to larger bodies balanced against gas drag: In this case we have: 

1/tR = 1/tgd: (44) 

where we use equation (37) for tR. This gives 



M2 



In (2 + Vs) - 1 



(45) 



Numerically for In(A^) = 1: 



\Cd^)\ j \b.2A\]) V r ) \Wm^Mi ) ' ^ ' 

Thus, even when mp has grown to the order of an Earth mass, the dispersion velocity ax 
plausibly remains such that u^.'^H/r, for r in the AU range, justifying the application of 
gas drag and maintaining a strong degree of gravitational focusing. For example, when 
m« = 10^® g, m^j = 1 Me and = 2.5 x IQ-^CdE, we find H^/r ~ 0.025, which gives 
an enhancement factor due to gravitational focusing (see above) of ~ 400 at 5.2 AU. 
This order of enhancement is needed in order to accumulate cores in the Earth mass 
range. 

1/5 

In addition, the velocity dispersion ax increases as and k, defined in the paragraph 
on runaway accretion, is equal to 14/15. Since this is less than unity, the accumulation 
does not occur in a runaway manner in which a single core outgrows everything else. 
Instead, a number of cores of comparable mass, or oligarchs, that are adequately sep- 
arated, grow together. But, nonetheless, gravitational focusing can and must remain 
important if core accumulation to the Earth mass range is to occur within the gas disk 
lifetime. 



— Scattering due to isolated larger bodies balanced against disk tides when gas is present: 
When all or most of the small mass planetesimals are accumulated, the likely situation 
that develops is one for which large mass cores, or oligarchs, of comparable mass mp are 
in near circular orbits separated by 10/ Hill radii, in isolation (Chambers et al 1996). 
The isolation mass may be estimated from the mean surface density of the accumulated 
solid material as: 



20/7rr% 
(3Me)V3 
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giving objects characteristically in the Earth mass range for a minimum mass solar 
nebula at 5.2 AU. But note that, for such a nebula, the isolation mass oc r^/^, being 

Such masses may slowly interact building up velocity dispersion and resulting in further 
accumulation (Papaloizou & Larwood 2000). For such large planetesimal masses mp, 
tidal interaction with the disk, should there be gas present, becomes more important 
than gas drag. Then we should use the timescale tg as given by equation (78) (see below) 
rather than tgd in equation (44) which determines the velocity dispersion equilibrium. 
The timescale tg is obtained from the protoplanet-disk interaction theory described 
below. As this timescale is significantly shorter than the associated possible type I 
migration timescale tj, we expect that a quasi-equilibrium is set up before significant 
migration occurs and that this then drives the evolution of the system, which will proceed 
under conditions of quasi-equilibrium. We recall the relaxation timescale for a central 
point solar mass potential in the form: 

_ 5 M| fH^Yf r \3/2 

- ln(A) TTE^r^m^ [~ ) [iMJ ) ^^^^ 

where A = (SM^/to/j) (i^/j/r)^ . Equating te from equation (78) and tf>, assuming 
e < H/r and = 1, gives the equilibrium semithickness of the protoplanet distribution 
in terms of the gas disk semithickness H: 



H ' ' V4E; VI AU. 

Since Hp/ H depends only very weakly on A, it is essentially independent of the proto- 
planet mass and it depends only weakly on other parameters. For the representative 
values mp = O.IM®, S^/(4E) = lO"-'^ and Hp/r = 0.01, we obtain Hp = 0.15H at 1 AU. 
Hence the protoplanet swarm is expected to remain thin and confined within the gaseous 
nebula as is confirmed by numerical simulations (Papaloizou & Larwood 2000). Note 
too that this implies that gravitational focusing will be important, enabling core mass 
build up to the Earth mass regime. However, disk tide induced type I orbital migration 
then becomes important on a 10^ yr timescale (Ward 1997, Tanaka et al 2002), threat- 
ening the survival of the cores against faUing into the central star (see below) . 



— The gas free case: When gas is absent so is the problem of inward core migration due 
to type I migration. However, should all the mass of solids be within the massive cores, 
there is nothing to balance velocity dispersion growth other than the effects of inter-core 
collisions. These are only likely to be limited when the relative velocity is comparable 
to the escape velocity (Safronov 1969). If the final phase of the accumulation of the 
terrestrial planets occurs in this way (Chambers Sz Wetherill 1998), we may estimate 
the evolution timescale for attaining 1 at 1 AU using equation (15), as gravitational 
focusing is of marginal importance, to be ~ 2 x 10^ yr. This timescale being oc implies 
an unacceptably long time to form Uranus and Neptune for minimum mass solar nebula 
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models. In addition, the escape velocity from a 1 M0 core exceeds the escape velocity 
from the sun if r > 12 AU, implying that escapers will be generated in this region. 
Goldreich et al (2004) suggest that a residue of small particles derived from collisions 
amongst the original planetesimals, such as has been produced in calculations directed 
towards explaining the distribution of Kuiper belt objects, by Kenyon & Bromley (2004), 
may remain and damp the velocity dispersions, reducing both evolution times and 
the orbital eccentricities of the remaining growing cores, once their number has been 
adequately reduced. 

3.2. Formation of giant planets: the core accretion model 

In order to produce a Jovian mass giant planet, the solid cores discussed above must 
accrete gas. According to the critical core mass model first developed by Perri and 
Cameron (1974) and Mizuno (1980; see also Stevenson 1982, Wuchterl 1995, Papaloizou 
& Terquem 1999), the solid core grows in mass along with the atmosphere in quasi- 
static and thermal equilibrium until it reaches the critical mass. Merit, above which 
no equilibrium solution can be found for the atmosphere. As long as the core mass is 
smaller than Merit , the energy radiated from the envelope into the surrounding nebula 
is compensated for by the gravitational energy which the planetesimals entering the 
atmosphere release when they collide with the surface of the core. During this phase of 
the evolution, both the core and the atmosphere grow in mass relatively slowly. By the 
time the core mass reaches Merit, the atmosphere has grown massive enough so that its 
energy losses can no longer be compensated for by the accretion of planetesimals alone. 
At that point, the envelope has to contract gravitationally to supply more energy. This 
is a runaway process, leading to the very rapid accretion of gas onto the protoplanet and 
to the formation of giant planets such as Jupiter. In earlier studies, it was assumed that 
this rapid evolution was a dynamical collapse, hence the designation 'core instability' for 
this model. The critical mass can be estimated to be in the range 5-15 depending 
on physical conditions and assumptions about grain opacity etc. (see, e.g.. Pollack et 
al 1996, Ikoma et al 2000, Papaloizou & Nelson 2005). There is thus a need to evolve 
beyond a typical isolation mass appropriate to a minimum solar nebula model and some 
authors have suggested that a nebula up to ten times more massive than the minimum 
mass solar nebula is needed (e.g.. Pollack et al 1996). 

However, this late evolutionary phase is not yet well studied because it depends on 
uncertain global phenomena such as inhibition of planetesimal accretion due to isolation, 
competition with other cores, gap formation (Thommes et al 2002) and how mobile the 
cores are under a modified type I migration which may be stochastic in the presence of 
turbulence (Nelson Sz Papaloizou 2004) or be reversed by large scale toroidal magnetic 
fields (Terquem 2003). Indeed, it is beheved that in order for giant planet formation 
to occur and to be able to fit the observed distribution of extrasolar planets, type I 
migration rates obtained for a laminar disk have to be reduced by such processes (e.g., 
Ahbert et al 2004). However, this has yet to be explored fully. 
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Time-dependent numerical calculations of protoplanetary evolution by Boden- 
heimer & Pollack (1986) and Pollack et al (1996) support the critical core mass model, 
although they show that the core mass beyond which runaway gas accretion occurs, 
which is referred to as the 'crossover mass', is slightly larger than Merit-, and that the 
very rapid gravitational contraction of the envelope is not a dynamical collapse. The 
designation 'crossover mass' comes from the fact that rapid contraction of the atmo- 
sphere occurs when the mass of the atmosphere is comparable to that of the core. Once 
the crossover mass is reached, the core no longer grows significantly through accretion 
of planetesimals. 

The evolutionary calculations by Pollack et al (1996) show that the evolution of 
a protoplanet with a gaseous envelope is governed by three distinct phases. During 
phase 1, runaway planetesimal accretion, as described above, occurs, bringing about 
the depletion of the feeding zone of the protoplanet leading to isolation. At this point, 
when phase 2 begins, the atmosphere is massive enough that the location of its outer 
boundary is determined by both the mass of gas and the mass of the core. As more gas 
is accreted, this outer radius moves out, so that more planetesimals can be captured 
enabling yet more gas to enter the atmosphere. The protoplanet grows in this way until 
the core reaches the crossover mass, which should be identified with the Mizuno (1980) 
critical core mass, at which point runaway gas accretion occurs and phase 3 begins. 
The timescale for planet formation is determined almost entirely by phase 2, i.e. by 
the rate at which the core can accrete planetesimals, and is found to be a few million 
years at 5 AU, being comparable to the gas disk lifetime. Note that although phase 3 is 
relatively rapid compared to phase 2 for cores as large as 15 , it may become much 
longer for smaller mass cores which become critical for reduced planetesimal accretion 
rates (see Pollack et al 1996, Papaloizou & Terquem 1999 and below). 

Note that isolation of the protoplanetary core, as it occurs at the end of phase 1, 
may be prevented under some circumstances by tidal interaction with the surrounding 
gaseous disk (Ward & Hahn 1995) or by interactions with other cores (Chambers et al 
1996). 

Conditions appropriate to Jupiter's present orbital radius arc normally considered 
and then for standard 'interstellar' dust opacities the critical or crossover mass is found 
to be around 15 M®. This is rather large as recent observations indicate that Jupiter 
has a sohd core not exceeding ~ 5 M®. 

The similarity between the critical and crossover masses is due to the fact that, 
although there is some liberation of gravitational energy as the atmosphere grows in 
mass together with the core, the effect is small as long as the atmospheric mass is small 
compared to that of the core. Consequently, the hydrostatic and thermal equilibrium 
approximation for the atmosphere is a good one for core masses smaller than the critical 
value. 

3.2.1. Equations governing an envelope at equilibrium: 



Planet formation and migration 



32 



We denote R the spherical polar radius in a frame with origin at the center of the 
protoplanet. Wc assume that we can model the protoplanct as a spherically symmetric 
nonrotating object. We consider the envelope at hydrostatic and thermal equilibrium. 
The equation of hydrostatic equilibrium is: 

where P is the pressure, p is the mass density per unit volume and g = GM{R)/R^ is 
the acceleration due to gravity, M{R) being the mass contained in the sphere of radius 
R (this includes the core mass if R is larger than the core radius). We also have the 
definition of density: 

-=4rfp. (51) 

At the base of a protoplanetary envelope, the densities are so high that the gas 
cannot be considered to be ideal. Thus we adopt the equation of state for a hydrogen 
and helium mixture given by Chabrier et al (1992). We take the mass fractions of 
hydrogen and helium to be 0.7 and 0.28, respectively. We also use the standard equation 
of radiative transport in the form: 

dT ^ -3kp L 

dR IGaTUnR^' ^ ' 

Here T is the temperature, k is the opacity, which in general depends on both p and 
T, (7 is the Stefan-Boltzmann constant and L is the radiative luminosity. Denoting the 
radiative and adiabatic temperature gradients by V^ad and Vad, respectively, we have: 

d\nT\ 3kL 

(53) 
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rad 
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(54) 



with the subscript s denoting evaluation at constant entropy. In writing equation (53), 
we have assumed that the only energy source comes from the accretion of planetesimals 
onto the core. The gravitational energy of the planetesimals going through the planet 
atmosphere is released near the surface of the core which, as a result, outputs a total 
luminosity Lcore given by (sec, e.g., Mizuno 1980; Bodenheimer & Pollack 1986): 

^core — ■ (,t)Oj 

^core 

Here Mcore and Vcore are respectively the mass and the radius of the core, and Mcore 
is the planetesimal accretion rate onto the core (note that some of the planetesimals 
entering the atmosphere may not end up colliding with the core). 

Note however that, as the the interiors of these protoplanet models are convective 
and approximately isentropic, the above formalism can be used when the energy 
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source corresponds to settling of the accreting gaseous atmosphere. One can use the 
conservation of mass and energy to determine the evolution of the protoplanet (see 
Papaloizou & Nelson 2005, and sections 3.2.3 and 3.2.4 below). 

When Vrad < Varf, the atmosphere is stable to convection and thus all the energy 
is transported by radiation, i.e. L = Lcore- When Vrad > Vad, there is instability to 
convection. Then, part of the energy is transported by convection, and Lcore = L+L^onv, 
where Lconv is the luminosity associated with convection. We use mixing length theory 
to evaluate Lconv (Cox & Giuli 1968). This gives: 



Lconv — T^Li CpAj^i 




1 3/2 
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J ^ 


^P9 





(56) 



where A^i = \amiP / {dP / dR)\ is the mixing length, ami being a constant of order unity, 
{dT/dR)g = V adT {dlnP/dR), and the subscript P means that the derivative has to be 
evaluated for a constant pressure. All the required thermodynamic parameters are given 
by Chabrier et al (1992). In the numerical calculations presented below we take ami = 1- 

— Inner Boundary: We suppose that the planet core has a uniform mass density 
Pcore- The composition of the planetesimals and the high temperatures and pressures 
at the surface of the core suggest pcore = 3.2 g cm~^ (e.g., Bodenheimer & Pollack 1986 
and Pollack et al 1996), which is the value we will adopt throughout. The core radius, 
which is the inner boundary of the atmosphere, is then given by: 

( 2,Mcore \"^ 

"^core = T— ■ (57) 

At i? = Tcore, the total mass is equal to Mcore- 

— Outer Boundary: In the models that will be discussed below, the outer boundary 
of the atmosphere will be taken to be either free or at the Roche lobe radius tl of the 
protoplanet given by: 

^^=oU^ ^' (58) 



3 V3M^ 

where Mpi — Mcore + Matm IS the planet mass, Matm being the mass of the atmosphere, 
and r is the orbital radius of the protoplanet in the disk. 



3.2.2. Models enveloped by the disk with luminosity derived from planetesimal accretion: 

These models must match to the state variables of the accretion disk at the sur- 
face of the Roche lobe. We denote the disk midplane temperature, pressure and mass 
density at the distance r from the central star by Tm, Pm and Pm, respectively. To 
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obtain the needed matching, we require that at i? = r^, the mass is equal to Mpi, the 
pressure is equal to and the temperature is given by: 



where we approximate the additional optical depth above the protoplanet atmosphere, 
through which radiation passes, by: 

Tl = (Pm, T^) PrnTL- (60) 

For a particular disk model, at a chosen radius r, for a given core mass Mcore and 
planctcsimal accretion rate Mcore, Papaloizou & Tcrqucm (1999) solve equations (50), 
(51) and (52) with the boundary conditions described above. The opacity, taken from 
Bell & Lin (1994), has contributions from dust grains assumed to be 'interstellar', 
molecules, atoms and ions. 

For a fixed Mcore at a given radius, there is a critical core mass Merit above which 
no solution can be found, i.e. there can be no atmosphere in hydrostatic and thermal 
equilibrium confined between the radii rcore and around cores with mass larger than 
Merit, as explained above. For masses below M^rit, there are (at least) two solutions, 
corresponding to a low-mass and a high-mass envelope, respectively. 

To estimate how large Mcore should be before a significant atmosphere can start to 
form, we suppose it to be isothermal, which is the most favourable situation for retaining 
a massive atmosphere. When the atmosphere starts to form, we may neglect its mass, 
so the gas density will be given by: 

/ GMcore \ 

We estimate that the atmosphere can start to be significant when 47rr j^^^gPfe/S = fMcore, 
where / < 1 is an estimated mass fraction contained in the atmosphere. This gives: 

/^2q2\3/2/ ^ \l/2/ /ofl7^ \\3/2 



1/4 





M_=^ i,^] (ln(=:i^)) . (62) 



Taking E = 200 g cm^, H/r = 0.05, r = 5.2 AU and / = 0.01, we obtain 
Meore = 0.064 Me. 

Figure 6 shows the total protoplanet mass Mpi as a function of core mass Mcore 
at different radii in a steady state disk with a = 10~^ (Shakura & Sunyaev 1973) and 
M = 10~^ Mq yr~^ (see Papaloizou & Tcrqucm 1999). In each frame, the different 
curves correspond to planetesimal accretion rates in the range IQ-i^-lQ-^ M® yr"^ 
The critical core mass is attained at the point where the curves start to loop backwards. 
It was demonstrated by Stevenson (1982) that the qualitative behaviour of the curves 
shown in Figure 6 can be produced using a simplified model with only radiative energy 
transport, ideal gas equation of state and constant opacity. It thus appears to result 
fom the core plus envelope nature of the models rather than details of the physics. 
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T„=2519.74 K, P„=29466.3 cgs T„=1968.16748 K, P„=18276.2773 cgs 
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Figure 6. Plots of total mass Mpi, in units M^, as a function of core mass Mcore, 
in units Mq, at different locations r in a steady state disk model with a = 10~^ 
and gas accretion rate M = 10~^ yr~^. From left to right and top to bottom, 
the frames correspond to r = 0.05, 0.06, 0.15, 0.5, 1 and 5 AU, respectively. The 
midplanc temperature and pressure at these locations are indicated above each frame. 
Each frame contains six curves which, moving from left to right, correspond to core 
luminosities derived from planetesimal accretion rates of Mcore = 10~^^, 10"^'^, 10~^, 
10~^, 10"'' and 10~^ yr~^, respectively. The critical core mass is attained when 
the curves first begin to loop backwards when moving from left to right. These plots 
are taken from Papaloizou & Terquem (1999). 
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When the core first begins to gravitationally bind some gas, the protoplanet is on 
the left on the lower branch of these curves. Assuming Mcore to be constant, as the 
core and the atmosphere grow in mass, the protoplanet moves along the lower branch 
up to the right, until the core reaches Merit- At that point, the hydrostatic and thermal 
equilibrium approximation can no longer be used for the atmosphere, which begins to 
undergo very rapid contraction. Figure 6 indicates that when the core mass reaches 
Merit, the mass of the atmosphere is comparable to that of the core, in agreement 
with Bodenheimer & Pollack (1986). Since the atmosphere in complete equilibrium is 
supported by the energy released by the planetesimals accreted onto the protoplanet, 
the critical core mass decreases as Mcore is reduced. 

For a = 10'^ and M = 10~^ Mq yr~^, the critical core mass at 5 AU varies 
between 16.2 and 1 as the planetesimal accretion rate varies between the largest 
and smallest value. The former result is in good agreement with that of Bodenheimer 
& Pollack (1986). Note that there is a tendency for the critical core masses to increase 
as the radial location moves inwards, the effect being most marked at small radii. At 
1 AU, the critical mass varies from 17.5 to 3 as the accretion rate varies between 
the largest and smallest value, while at 0.05 AU these values increase still further to 42 
and 9 M®, respectively. 

In Figure 7 wc plot the critical core mass M^^n versus the location r for three 
different steady disk models. These models have a = 10^^ and M = 10^'' M0 yr~^, 
a = 10^2 ^ ;lq-8 y^-i^ ^^^i ^ ^ iQ-^ ^ iQ-a y^-i^ respectively 

Here again, in each frame, the different curves correspond to planetesimal accretion rates 
in the range lO^^^-lO^® yr^^. Similar qualitative behavior is found for the three disk 
models, but the critical core masses are smaller for the models with M — 10~^ M0 yr~^, 
being reduced to 27 and 6 at 0.05 AU for the highest and lowest accretion rate, 
respectively. 

3.2.3. Models enveloped by the disk with no planetesimal accretion: 

For these models, subsequently denoted as of type A, just as for those with plane- 
tesimal accretion as energy source, wc assume the structure extends to the Roche lobe 
or boundary of the Hill sphere beyond which material must be gravitationally unbound 
from the protoplanet. There it matches to the external disk. Thus the boundary con- 
ditions for these models are the same as those with planetesimal accretion as energy 
source. 

Papaloizou & Nelson (2005) adopt disk parameters appropriate to 5 AU from the disk 

model of Papaloizou & Tcrquem (1999) with a = lO^"' (Shakura and Sunyaev 1973) 
and steady sate accretion rate of 10"'' yr~^. Accordingly, — 140.047 K and 
Pm = 0.131 dyn cm^^. 

For models of type A, specification of the core mass and the total mass Mpi is enough 
to enable it to be constructed. This is because, although the luminosity is apparently 
an additional parameter, the degree of freedom associated with it is specified by the 
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Figure 7. Critical core mass Merit, in units Mq, vs. location r, in units AU, in 
a steady state disk model with viscous parameter a and gas accretion rate M. The 
different plots correspond to a = 10~^ and M = 10"'' Mq yr~^ (a), a = 10~^ and 
M = 10-* Mq yr-i (b), and a = 10"^ and M = 10"* Mq yr"! (c). Each plot 
contains six curves which, moving from bottom to top, correspond to core luminosities 
derived from planetesimal accretion rates of Mcore = 10~^^, 10~^°, 10~^, 10~^, 10~^ 
and 10~^ yr~^, respectively. These plots are taken from Papaloizou & Terquem 
(1999). 
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requirement that the radius of the planet atmosphere should be the Hill radius which, 
for a fixed orbital radius, depends only on Mpi. Thus, the luminosity and the total 
energy of the protoplanet may be written as L = L{Mpi) and E = E[Mpi). 
For small changes in mass, the change in E is dE = {dE / dMpi)dMpi. If the change 
occurs over a time interval dt, we must have: 

dMpi dt ^ ' 

This determines how the mass and hence the luminosity of the protoplanet vary with 
time and therefore the evolutionary track is specified. 

3.2.]^.. Models with free boundary accreting from the protostellar disk: 

In contrast to embedded models of type A, we can consider models that have bounda- 
ries detached from and interior to the Roche lobe which still accrete material from the 
external protoplanetary disk that orbits the central star. Detachment from the Roche 
lobe and the development of a free boundary are expected once the planet goes into 
a rapid accretion phase. This is also expected because numerical simulations of disk- 
planet interactions have shown that, once it becomes massive enough, a protoplanet 
forms a gap in the disk but is still able to accrete from it through a circumplanetary 
disk (see, e.g., Bryden et al 1999, Kley 1999, Lubow et al 1999, R. Nelson et al 2000). 
We thus consider models with free boundaries which are able to increase their mass 
and liberate gravitational energy through its settling. These models are referred to as 
type B. 

In contrast to models of type A, for an assumed externally supplied accretion rate, 
models of type B form a two parameter family in that, without specification of the 
energy source, and given their freedom to determine their own radius, they require 
specification of both Mpi and L in order for a model to be constructed. Thus 
E — E{Mpi,L). Accordingly, for small changes in mass and luminosity, the change 
in E is dE^ {dE / dMpi)dMpi + {dE/dL)dL. 

Now, for these models, matter is presumed to join the protoplanet on its equator after 

having accreted through a circumplanetary disk. In this case, we assume the accretion 

rate to be prescribed by the dynamics of the disk-planet interaction while gap formation 

is taking place. This is in agreement with simulations of disk-planet interactions where 

it is found that an amount of material comparable to that flowing through the disk 

may be supplied to the protoplanet (Bryden et al 1999, Kley 1999, Lubow et al 1999, 

R. Nelson et al 2000 and simulations presented in section 4). In arriving at the planet 

equator, all available gravitational binding energy of —GMpi/rg per unit mass, being 

the surface radius, has been liberated and so an amount of energy —GMpidMpi/rg must 

be subtracted from dE in order to obtain the energy available to replace radiation losses. 

Therefore, if the changes occur over an interval dt^ we must have: 

GM^dM^ dE dE GM^dM^ 

dE + = ^dMp, + -dL + = -Ldt. (64) 
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Thus, total energy conservation for models of type B enables the calculation of 
evolutionary tracks through: 



dE ^ GMpi 



pi 



dMpi dEdL 



+ ^-J7 = -L- (65) 



dt dL dt 



Note that, as we regard the accretion rate dMpi/dt as specified for these models, 
equation (65) enables the evolution of L to be calculated. 

Thus, equations (63) and (65) constitute the basic equations governing the evolution of 
models of type A and type B, respectively. 

Note that any input from planctcsimal accretion during and after the phase when 
the core becomes critical is neglected. Wc recall that after this point planctcsimal 
accumulation proceeds by oligarchic growth (Ida & Makino 1993, Kokubo & Ida 1998). 
N-body simulations of protoplanetary core formation indicate that obtaining cores of 
a high enough mass ~ 15 to obtain evolutionary timescales within the gas disk 
lifetime is not an easy task to achieve during the oligarchic growth phase. This is in 
part due to planetary cores of a few repelling the surrounding planetesimals and 
opening gaps in the planctcsimal disk (e.g., Thommes et al 2002). In addition, it is 
reasonable to suppose that the formation of a series of neighboring cores may deplete 
planetesimals in the local neighborhood. We also note that if the slow down in evolution 
rate due to the long thermal timcscalc of the atmosphere can be obviated, a depletion 
of planetesimals brings about a favourable situation for shortening the evolutionary 
timescale as a potential energy source for sustaining the models without gas accretion is 
removed. The thermal timescale may be reduced by reducing the contribution of grain 
opacity from the 'interstellar value' (see figure 8, Ikoma et al 2000 and Papaloizou & 
Nelson 2005). If the grain opacity in the outer layers of the atmosphere can be reduced 
by a factor of 100, the evolutionary timescale for a protoplanet with a 5 core can 
be reduced from 3 x 10^ yr to a value within the expected gas disk lifetime. Such a 
reduction may occur either because of grain depletion or because the grains have been 
accumulated into particles of larger sizes than the original (Ikoma et al 2000). 

However, regardless of opacity, it is found that the evolutionary timescales for the 
type A models decrease rapidly once the gas mass is about twice the core mass (see 
figure 8). At this point, detachment from the Roche lobe and transition to a model of 
type B is expected (Papaloizou & Nelson 2005). Some of these are illustrated in figure 9. 
It is found that these models rapidly attain a radius <2 x 10^*^ cm for any reasonable 
accretion rate from the disk. Furthermore, these protoplanets may attain luminosities 
of - 10-3-10"^ L© for a time - 10^-10^ yr. 



3. 3. Form,ation of giant planets: gravitational instabilities 

Gravitational instabilities in protostellar disks have been considered as a mechanism for 
forming giant planets since Kuiper (1949) proposed that planets were formed through 
gravitational breakup of gaseous rings. The advantage of this mechanism is that it leads 
to giant planets on a very short timescale, on the order of the dynamical timescale. This 



Planet formation and migration 



40 




Figure 8. This figure illustrates the evolution of protoplanet models of type A 
which maintain contact with the protoplanetary disk and fill their Roche lobes while 
they accrete from it. They have fixed solid core masses of 15 and are situated 
at 5 AU. The upper left panel shows the luminosity, in cgs units, as a function of 
their increasing mass M, in cgs units. The upper right panel shows the gas accretion 
rate, in y""'^, as a fimction of mass, while the lower left panel gives the accretion 
time M/M, in yr, as a function of mass. The lower right panel gives the temperature 
Teffp, being the radiation temperature on the Roche lobe required account for the 
luminosity of the protoplanet, as a function of mass. The model illustrated with 
the full curve is embedded in a standard disk while the model illustrated with a 
dotted curve is embedded in a disk with the same temperature but with a density 
ten times larger. These two protoplanet models with standard opacity show very 
similar behavior indicating lack of sensitivity to the detailed boundary conditions. In 
addition, we illustrate two models with this core mass embedded in a standard disk 
but with opacities which have a reduction factors of ten and one hundred (dashed 
curves and triple dot dashed curves, respectively) that are constant for T < 1600 K 
and which then decrease linearly to unity at T = 1700 K. Note the rapid evolutionary 
speed up that occurs when these models attain ~ 35 Mq. At this point, detachment 
from the Roche lobe and transition to a model of type B is expected. (These plots are 
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Figure 9. This figure illustrates the evolution of protoplanet models of type B which 

accrete from the protoplanctary disk at an assumed rate of one Jupiter mass in 9x 10'* yr 
but arc detached from their Roche lobes. They have fixed solid core masses of 15 
and are situated at 5 AU. The upper left panel shows the planet mass, in cgs units, as 
a function of time, in yr. The upper right panel shows the protoplanet radius, in cgs 
units, as a function of time. The lower left panel gives the luminosity, in cgs units, as a 
function of time. The lower right panel gives the effective temperature as a function of 
time. The four models shown correspond to differing initial conditions corresponding 
to different starting masses and luminosities. The resulting evolutionary tracks tend 
to show convergence as time progresses. (These plots are taken from Papaloizou & 
Nelson 2005). 
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is an important argument, as we have seen above that the build-up of a massive core 
may be prohibitively long in some parts of the protostellar nebula. Also, some of the 
extrasolar planets discovered so far are very massive, with a mass several times that of 
Jupiter. It is not clear that these objects can form through the core-accretion scenario, 
as gap formation inhibits gas accretion above some critical mass (see below). 

It was initially beheved that axisymmetric instabihties in disks would lead to the 
formation of rings that would further break-up into condensations. Using the theory of 
viscous accretion disks developed by Lynden-Bell and Pringle (1974), Cameron (1978) 
calculated evolution models of the primitive solar nebula which were unstable against 
axisymmetric perturbations. 

The parameter governing the importance of disk self-gravity is the Toomre 
Q parameter which, in a Keplerian disk, can be written as Q ~ i7cs/(7rG'S) ~ 
M^H/{Mdr), with being the disk mass contained within radius r, H the disc 
semithickness and Cg the sound speed. Typically H/r ~ 0.1 for protostellar disks 
(Stapelfeldt et al 1998) such that the condition for the importance of self-gravity, 
(5 ~ 1, gives Mrf ~ O.lMy,. The characteristic scale associated with growing density 
perturbations in a disk undergoing gravitational instability with Q ~ 1 is on the order 
of H, and the corresponding mass scale is Md{H/rY ~ Mi,{H/rY, which is on the order 
of one Jupiter mass for H/r ~ 0.1 and = 1 M©. 

There are, however, problems with forming planets through gravitational ins- 
tabilities. Fragmentation, which follows the development of axisymmetric modes, 
requires Q < 1, whereas non axisymmetric modes develop for values of Q slightly larger 
than unity. As noted by Laughlin & Bodenheimer (1994), it is unlikely that a disk would 
ever get massive enough for Q decrease below unity. The disk builds up in mass as a 
result of the collapse of an envelope. The Q parameter therefore decreases with time from 
high values. When it reaches values slightly larger than unity, non axisymmetric (spiral) 
modes develop. These lead to outward angular momentum transport (Papaloizou & 
Savonije 1991; Heemskerk et al 1992; Laughlin & Bodenheimer 1994; Pickett et al 1998) 
that results in additional mass growth of the central star. Therefore Q is decreased and 
gravitational stability is restored. Since the redistribution of mass may occur on the 
dynamical timcscalc (a few orbits) of the outer part of the disk, it follows that Q can 
never decrease to below unity unless the disk can be cooled down (so that decreases) 
or matter can be added (so that E increases) on the dynamical timescale, which is 
unlikely to happen in the regions of planet formation (it could happen however in the 
optically thin parts of the disk, typically beyond ~ 40 AU) . 

Nonetheless, condensations may form even if the disk is unstable only against non- 
axisymmetric perturbations, providing the density gets high enough in the spiral arms 
(Boss 1998, Mayer et al 2004). There are indications, however, that the clumps that 
appear in the numerical simulations of gravitationally unstable disks are not long- 
lived (see, e.g., Durisen 2001 and references therein). The outcome of gravitational 
instabilities is very sensitive to the thermal physics of the gas. Chemical composition 
dependence is expected but the correlation of the existence of planets with metallicity 
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has yet to be demonstrated. When the equation of state is adiabatic, heating of the gas 
prevents the formation of clumps. In contrast, condensations form in an isothermal disk. 
But if the outer radius of the disk is allowed to expand under the effect of gravitational 
torques, these clumps tend to disappear, dissolving back into the gas (Pickett et al 
1998, 2000a, b; Boss 2000). Simulations including heating via dynamical processes and 
radiative cooling produce spiral structures that are less distinct than in isothermal disks 
and that do not collapse into condensed objects (A. Nelson et al 2000). Although the 
question of whether the clumps can survive or not is still being debated, it seems unlikely 
that Jupiter mass planets can be formed through gravitational instabilities in optically 
thick regions of disks. However, as noted in section 3.1.2, the accumulation of solid 
particles into planetesimals may be assisted by gravitational instabilities in the gas disk 
through particle accumulation at pressure maxima in spiral structures (Rice et al 2004, 
Durisen et al 2005). 

Note that fragmentation may also occur before a disk is completely formed, during 
the initial collapse of the protostellar envelope. Such opacity limited fragmentation has 
been estimated to produce objects with a lower mass limit of 7 Jupiter masses (Low & 
Lynden-Bell 1976, Masunaga & Inutsuka 1999), but there is no definitive argument to 
rule out somewhat smaller masses (Bodenheimer et al 2000). It is possible that both a 
disk and fragments may form simultaneously out of the protostellar envelope, the relative 
importance of the two processes depending for instance on the angular momentum 
content of the envelope, on the strength of any magnetic field (so far neglected in 
disk fragmentation calculations) and possibly on the initial dumpiness. Large scale 
observations of class envelopes so far do not rule out the presence of clumps with 
masses smaller than about 10 Jupiter masses (Motte & Andre 2001). To explain the 
large eccentricities of many of the extrasolar planets, Papaloizou & Terquem (2001) 
and Terquem & Papaloizou (2002) have considered a scenario in which a population of 
planetary mass objects are assumed to form rapidly through a fragmentation process 
occurring in a disk or protostellar envelope on a scale of 100 AU (see section 5 below). 

4. Disk— planet interactions and orbital migration 

The discovery of the 'hot Jupiters' led to the realization of the importance of orbital 
migration. This comes about because of the perceived difficulties in accumulating a core 
directly at small radii interior to the so called 'snow line' (Bodenheimer et al 2000). 
Instead, it is possible the protoplanet formed at larger radii and migrated inwards. 
So far, three mechanisms have been proposed to explain the location of planets at very 
short orbital distances. One of them relies on the gravitational interaction between two 
or more Jupiter mass planets, which may lead to orbit crossing and to the ejection of 
one planet while the other is left in a smaller orbit (Rasio & Ford 1996, Weidenschilling 
& Marzari 1996). However, this mechanism cannot account for the relatively large 
number (about 20%) of short-period planets observed. Another mechanism is the so- 
called 'migration instabihty' (Murray et al 1998, Malhotra 1993). It involves resonant 
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interactions between the planet and planetesimals located inside its orbit which lead to 
the ejection of a fraction of them while simultaneously causing the planet to migrate 
inwards. To move a Jupiter mass planet from about 5 AU down to very small radii 
through this process, a disk containing about 1 Jupiter mass of planetesimals, and thus 
about 0.1 Mq of gas, inside 5 AU is required. Such a massive disk is unlikely and 
furthermore it would be only marginally gravitationally stable. The third mechanism, 
that we are going to focus on here, involves the tidal interaction between the protoplanet 
and the gas in the surrounding protoplanetary nebula (Goldreich & Tremaine 1979, 1980; 
Lin & Papaloizou 1979, 1993 and references therein; Papaloizou & Lin 1984; Ward 1986, 
1997b). 

Below we review the possible main types of orbital migration associated with this tidal 
interaction. These are (i) type I migration, which applies to an embedded protoplanet of 
small mass for which the disk response can be considered using linear analysis (e.g., Ward 
1997), (ii) type II migration, which applies when the protoplanet is massive enough to 
open a gap (Lin & Papaloizou 1986) and (iii) runaway or type III migration (Masset Sz 
Papaloizou 2003, Artymowicz 2004), which is a new form of potentially fast migration, 
applicable to massive disks, that could be driven by coorbital torques. 
Orbital migration in each of these cases is induced by the angular momentum exchange 
with disk material that occurs through disk-planet interaction. This may be through 
wave excitation, as in type I migration, a combination of wave excitation and shock 
dissipation, as occurs in type II migration, or direct exchange with disk material 
traversing the orbit, as has been suggested for type III migration. 
In the appendix we give a simple description of the global angular momentum balance 
and conservation when a protoplanet in circular orbit is embedded in a gaseous disk 
with which it interacts gravitationally and consequently migrates radially. We show 
that the angular momentum exchange with the protoplanet can be accounted for by 
considering viscous and wave fluxes of angular momentum, that can be measured a long 
way from the protoplanet, as well as the difference between the rate of advection of 
angular momentum by gas flowing through the orbit and the rate of increase in the 
angular momentum content of regions in the neighborhood of the protoplanet because 
of non steady conditions. The angular momentum carried in wave/viscous fluxes is 
associated with type I and type II migration, whereas the advected flux is associated 
with the more recently proposed type III migration. 

We consider the timescales for the various migration processes and the potential 
problems introduced for planet formation theory. 

4.1. Angular momentum transfer due to scattering and type II migration 
4.1.1. The impulse approximation applied to a protoplanet in a gap: 

Here, we adopt a very simple picture of the angular momentum exchange process 
between a protoplanet in circular orbit and the disk. We work in a frame in which the 
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protoplanet appears at rest at the origin. For a protoplanet in circular orbit, this frame 
would have to be in a state of uniform rotation and we should include centrifugal and 
Coriolis forces to describe the motion of the disk matter. However, we shall neglect these 
and make a local approximation, supposing that a particle of disk matter approaches 
the protoplanet from large distance moving in a sraight line with impact parameter a 
which we shall take to be much less than the orbital radius of the protoplanet (Lin & 
Papaloizou 1979). 

On passing the protoplanet, the gravitational interaction produces an angle of deflection 
S between the asymptotic tanjents to the trajectory of the passing element of disk matter. 
In the local approximation, the disk matter recedes from the protoplanet along a straight 
line inclined to the line of approach by an angle S. 

We adopt local Cartesian coordinates with origin at the center of mass of the 
protoplanet. The x-axis points radially outwards, the y-axis points in the azimuthal 
direction in the direction of rotation while the ^-axis points in the vertical direction. 
We may assume that the initial trajectory is along x — a approaching from y — — oo 
with relative speed u. This corresponds to a disk with circular streamlines. 
As a result of the interaction, a velocity component Vx, parallel to the x-axis, is induced. 
This can be found from the x-component of the equation of motion in the form: 

where Mp is the planet mass. Assuming the interaction to be weak, we set x = a, 
y — ut and integrate with respect to time between — oo and +00. Hence, the final 
induced x-velocity is: 

r... r dt 2GMp 

Vx = -GMpaj_^ PTFFF ^ 

The angle of deflection is thus simply: 



6 



2GMp 



Associated with this interaction is a transfer of angular momentum between the disk 
matter and protoplanet. To see this, we first recall that, as the interaction between 
protoplanet and disk is conservative, energy conservation implies that the ultimate speed 
of recession is also u. Accordingly, there must have been a parallel momentum transfer 
per unit mass of magnitude A-u = u{l — cos 5) ~ u5'^/2. If the orbital radius of the 
protoplanet is R, the associated angular momentum transferred is: 

2RG^Ml , ^ 

A J = Ru5 = (66 

We should of course note that this transfer is directed along the original direction of 
motion as seen by the protoplanet. Thus, for disk matter on circular orbits interior to 
the protoplanet, the transfer is positive, while for exterior disk matter, the transfer is in 
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the opposite sense and negative. The interaction is accordingly frictional, with interior 
matter speeding up the protoplanet and exterior matter slowing it down. 
Note too that the trajectory prior to the local scattering interaction is assumed 
linear (corresponding to circular orbit when viewed globally). Subsequent returns of 
disk matter on circular orbit are required in order to make the frictional interaction 
persistent. The scattering itself, of course, disturbs that situation, so that there is an 
implicit assumption that dissipative or other processes work to restore circular orbits 
for returning trajectories. The indications from many numerical simulations are that 
disk viscosity is able to provide such an effect. 

The simple calculation of AJ given above neglects the fact that the calculation should 
be done in a rotating frame. However, this has been incorporated in calculations by 
Goldreich & Tremaine (1980) and Lin & Papaloizou (1993). As these calculations add 
nothing new with regard to the essential physics, we refer the reader to those works for 
the details. The end result is that the expression given in equation (66) is multiplied by 
a correction factor 

Co = -[2Xo(2/3) + Xi(2/3)]% 

where Kq and Ki denote modified Bessel functions. Thus: 

2RG^Ml , ^ 

AJ = Co^^. 67 

4- 1-2. Angular momentum exchange rate with the disk: 

We may use equation (67) to evaluate the total angular momentum exchange rate 
with the disk. To do this, we replace the relative speed u by u = R\Q — uj\, where u 
and Q are the angular rotation rates of the protoplanet and disk matter, respectively. 
Given that we have implicitly assumed the interaction occurs close to the protoplanet, 
we adopt the first order expansion |f2 — a;| = 3aa;/(2i?). 

Furthermore, each disk particle orbiting a distance a from the protoplanet suffers 
impulses separated by a time interval 27r/|Q — We can thus write down the rate 
of transfer of angular momentum between the disk matter interior to the protoplanet 
orbit and the protoplanet by integrating over the disk mass as: 

dJ_ 
H ~ 

This gives: 



KL/^J\VL - uj\da. (68) 



dt 27a;2 



3 



(69) 



Here, amin is the minimum distance between disk matter and the protoplanet orbit. 
Thus, there is an inherent assumption of a gap in the disk, or lack of coorbital material, 
which would have to be taken into account if the protoplanet is embedded in the disk and 
the interaction is linear. Therefore, the above analysis is applicable only in a nonlinear 



Planet formation and migration 



47 



regime when there is a significant gap, corresponding to type II migration. Note that 
the sign is positive for the more rapidly moving inner disk material. For the outer disk 
material, one can see from symmetry considerations that the same expression applies 
but with a sign reversal because it is more slowly moving than the protoplanet and thus 
attempts to slow it down. 



4.1.3. Gap formation: 



A gap forms because angular momentum exchange occurs locally in the disk, on each 
side of the protoplanet, in such a way as to cause the disk material to move away from 
the protoplanet. 

In order for a gap to be maintained, we should expect that the angular momentum 
exchange rate given by equation (69) be able to balance the angular momentum flow 
rate due to viscosity. The viscous flow rate is given by: 




= 37ri/Ei?^a;, (70) 

where u is the kinematic viscosity and E is the surface density of gas in the disk. Of 
course, we recall that u is not a real 'viscosity' but a way of representing the flow of 
angular momentum resulting from turbulent stresses. From what we have indicated 
above, we expect the condition for gap formation to be: 

(71) 




This is equivalent to: 

where is the mass of the central star, around which the planet orbits. To work with 
equation (72), we need to know what value of amin to use. Considerations of the pressure 
response in type I migration (Goldreich & Tremaine 1980, Artymowicz 1993) suggests 
that there should be a torque cut-off leading to a^m ~ H- This is also a statement 
that pressure effects are able to smooth out smaller scale perturbations in the vicinity 
of the protoplanet. A similar conclusion was reached by Papaloizou & Lin (1984) who 
considered the linear response of a protoplanet orbiting in a disk already with a gap 
with lengthscale H, which was found to be the physical scale that led to the maximum 
angular momentum exchange rate. 

However, the condition for strong nonlinearity, which is required for local gap formation, 
is that the radial lengthscale of the flow cannot be less than the Hill radius, r^, 
which being the largest possible scale of signigicant influence by the protoplanet 
gravitational field (Lin & Papaloizou 1979), must exceed H. Accordingly, we should have 
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Ornin > ^H- Taking ajnin — '^f'H and absorbing uncertainties in the constant coefficient 
Co, equation (72) gives: 



3MpCo 



> (73) 



8l7rM^ R^uj 
Taking Cq — 2/3, this leads to: 

^ > i^. (74) 

This condition is given by Lin & Papaloizou (1986, 1993). But note that it is only 
expected to be accurate to within a constant of order unity. For a typical disk model, 
v/^R^VL) 10~^, and we expect gap formation for planets slightly above the mass of 
Saturn. 

However, note that Rafikov (2002) argues that gap formation may occur for planets of 
mass as low as 2 M®. So far, gap formation at such low masses has not been seen in 
simulations (e.g., D'Angelo et al 2003). This may be because Rafikov (2002) suggested 
that a gap can be produced as the density waves excited by a planet become weak 
shock waves and dissipate. However, this occurs several scale heights away from the 
planet, beyond the region where a gap would be expected. Furthermore, the angular 
momentum deposition from the waves also extends over a region of several scale heights, 
being wider than an expected gap. In this context, note that an annulus of width 5H 
at 5.2 AU in a minimum mass solar nebula contains a Jupiter mass, so that a planet of 
several earth masses could not redistribute its angular momentum without significant 
migration. 

4.1.4- Type 11 migration: 

When the conditions for gap formation are satisfied, the outward flux of angular 
momentum due to viscous transport in the outer disk section is supplied by the planet. 
Similarly, the outward flux from the inner section is taken up by the planet. When 
these contributions do not balance, the planetary orbit changes its angular momentum 
and thus migrates. Normally, simulations indicate that the outer disk has the dominant 
effect, causing inward migration. This is illustrated in figure 10, taken from R. Nelson et 
al (2000), which shows the evolution of a planet of initially 1 Mj in a standard disk. In 
this case, the protoplanet was allowed to accrete as it migrated inwards on the viscous 
timescale of ~ 10^ yr. Inward migration on the viscous timescale R"^/ v is characteristic 
of type 11 migration driven by the evolution of the disk when the protoplanet mass 
does not exceed the local disk mass. In that case, using the standard a-prescription by 
Shakura and Sunyaev (1973), the migration timescale is given by: 

If we adopt Hjr = 0.1 and consider a in the range 10 '^-10 ^, we get t\i ~ 10^ or 10^ yr 
at i? = 1 or 5 AU, respectively. These timescales are much shorter than the disk lifetime 
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Figure 10. This figure shows the time evolution of the disk surface density contours 
and orbital radius as a function of time measured in orbital periods for a protoplanet 
with initial mass of 1 Mj. As the protoplanet undergoes type II migration inwards, 
its mass grows to 3.5 Mj under an assumption of maximal accretion. (The plots are 
taken from R. Nelson et al 2000). 



or estimated planetary formation timescales. The expression of tn is independent of the 
protoplanet's mass and the surface mass density of the disk, but it is imphcitly assumed 
here that the mass of gas within the characteristic orbital radius of the planet is at least 
comparable to the mass of the planet itself. If the disk is significantly less massive, the 
inertia of the planet slows down the migration (Syer & Clarke 1995, Ivanov et al 1999). 

4.. 1.5. Stopping type II migration: 

— Stopping type II migration at small radii: Planets subject to type 11 migration 
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would stop before falling onto the star if they entered a magnet ospheric cavity. Tidal 
interaction with a rapidly rotating star would also halts planet orbital decay at a few 
stellar radii (where the interaction becomes significant; see Lin et al 2000 and refer- 
ences therein). Both of these mechanisms have been put forth to account for the present 
location of the planet around 51 Pegasi (Lin et al 1996), and to explain the location of 
hot Jupiters more generally. 

A planet overflowing its Roche lobe and losing part of its mass to the central star would 
also halt at small radii. This is because during the transfer of mass the planet moves 
outward to conserve the angular momentum of the system (Trilling et al 1998). The 
planet stops at the location where its physical radius is equal to its Roche radius. Re- 
cent observations of atomic hydrogen absorption in the stellar Lyman a line during 
three transits of the planet HD209458b suggest that hydrogen atoms are escaping the 
planetary atmosphere (Vidal-Madjar et al 2003). 

— "The last of the Mohicans"...: The mechanisms reviewed above cannot account for 
the presence of giant planets orbiting their parent star at distances larger than about a 
tenth of an AU. 

It has been suggested that migration of a giant planet could be stopped at any radius 

if migration and disk dissipation were concurrent (Trilling et al 1998, 2002). If the 
disk dissipates while migration is taking place, then the drift timescale may increase in 
such a way that the planet stalls at some finite radius. Note however that this requires 
very fine tuning of the parameters (disk mass, disk lifetime etc.), as for a given disk 
mass the migration (viscous) timescale given by equation (75) decreases as the orbital 
radius decreases. Also, a major problem with this mechanism is to explain how the disk 
dissipates. Within this scenario, there is initially enough gas in the disk to push the 
planet down to some orbital radius. For typical disk parameters, only part of this gas 
may be accreted by the planet or leak through the gap to be accreted onto the star 
(Bryden et al 1999, Kley 1999, Lubow et al 1999, R. Nelson et al 2000). It is therefore 
not clear how the gas disappears. 

A giant planet could survive if after it formed there were not enough material left in 
the disk for significant migration to occur. It has been suggested that a scries of giant 
planets could actually assemble in the disk and disappear onto the star (Gonzalez 1997, 
Laughlin & Adams 1997). Then at some point the disk mass may be such that one more 
planet can be formed but not migrate (Lin 1997). This survivor is sometimes refereed 
to as the last of the Mohicans. 

Note that the prediction that the more massive giant planets that arc able to form gaps 
should be predominantly at larger radii is supported by the observations. 

— Planets locked in resonances: Orbital migration occuring for different planets at 
different rates, dependent on local disk parameters, can produce convergent migration 
and locking on to mean motion commensurabilities in a similar manner as is believed to 
occur as a result of tidally induced migration of sateUites in the solar system (Goldreich 
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1965). 

Pairs of extrasolar planets have been discovered in such commensurabihties. Ghese 876 
(Marcy et al 2001) and HD 82943 (Mayor et al 2004) are in 2:1 commensurabihties 
while 55 Cancri (McArthur et al 2004) exhibits a 3:1 commensurability. The origin of 
these commensurabihties is satisfactorily accounted for by disk planet interaction (e.g., 
Snellgrove et al 2001, Lee & Peale 2002, Papaloizou 2003, Kley et al 2004). However, 
there is no indication that the existence of these commensurabihties slowed down or 
halted inward migration. On the other hand, it is possible that for some planet masses 
and disk conditions, resonant trapping of planets could lead to a reversal of type II 
migration. This has been suggested to occur when a giant planet (e.g., Jupiter) migrates 
inwards and captures a lighter Saturn mass planet in 2:3 resonance (Masset & Snellgrove 
2001). However, it is unclear how often and for how long such situations can be realized. 
We also comment that apsidal resonances in which the apsidal lines of two planet 
orbits are aligned, but not commensurable, could also be produced through disk planet 
interactions which are dissipative and therefore tend to select a particular normal mode 
involving non zero eccentricity of the planetary and disk orbits jointly, while the latter 
is present (e.g., Papaloizou 2003). 

4-2. Linear torques and type I migration 

In this case, one is concerned with protoplanets that are not massive enough to form a 
gap in the disk. We are therefore interested in calculating the hnear response of the disk 
to an embedded orbiting protoplanet and using that to calculate the angular momentum 

exchange rate with the protoplanet. 

The perturbation exerted by the protoplanet on the disk propagates as density waves 
outside the Lindblad resonances and is evanescent inside these resonances, in the 
corotation region. The protoplanet exerts a torque on the density waves, which, 
together with the torque exerted at corotation, is responsible for the exchange of angular 
momentum between the disk rotation and the planet's orbital motion. Away from the 
Lindblad resonances, the waves have a small wavelength, so that the contribution to 
the net torque is small. Therefore, most of the torque is exerted in the vicinity of the 
Lindblad resonances, where the perturbation has a large wavelength, or in the vicinity 
of the corotation resonance (sec Tcrqucm et al 2000 and references therein). 
Thus, the torque exerted by the protoplanet on the disk can be calculated either by 
performing a full linear analysis of the excited waves in two or three dimensions (see, 
for example, Korycansky & Pollack 1993, Tanaka et al 2002, Tcrquem 2003, Tanaka 
& Ward 2004) or by summing up the contribution of point like torques exerted at 
Lindblad and corotation resonances (see, for example, Goldreich & Tremaine 1979, 
Artymowicz 1993, Ward 1997, Papaloizou & Larwood 2000). The results of fuU 
numerical simulations in general confirm the linear analysis (e.g., D'Angelo et al 2003, 
Bate et al 2003). Both analytical and numerical calculations find a positive (negative) 
torque on the disk (protoplanet) when the disk is laminar, circular and does not contain 
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any magnetic field. This corresponds to inward migration of tlie protoplanet. 
Tlie conditions for non linearity, or gap formation, have been explored by Papaloizou 
et al (2004) using both local shearing box simulations and fully global simulations. A 
sequence of local shearing box simulations, incorporating MHD turbulence with zero 
net magnetic flux, which have increasing MpE? / (M^H^), is illustrated in figure 11. The 
transition from being embedded to opening a gap is seen to occur when this parameter 
is around unity and it corresponds to the nonlinearity condition mentioned above that 
the radius of the Hill sphere should exceed the disk thickness. When this condition is 
satisfied, the viscous transport condition given by equation (74) is also satisfied in these 
simulations. 

4-2.1. Timescales: 

The resulting dependence of the linear (inward) migration time ti and eccentricity 
damping time te, calculated for a laminar inviscid disk as a sum of contributions from 
Lindblad resonances, on the gas disk thickness is such that tj = \J/j\ oc {H/RY and 
te = |e/e| oc {H/Ry. This scahng results in te being much shorter than tj. Accordingly, 
we expect circular orbits to be set up quickly and then migrate. But note that this 
may not be the case when the disk is eccentric. Then the disk interaction is expected 
to lead to a finite eccentricity in the protoplanet orbit (Papaloizou 2002). Furthermore, 
the migration direction may be affected. 

Tanaka et al (2002) performed a local three dimensional calculation for a 
planet in circular orbit embedded in an isothermal axisymmetric disk that summed 
contributions from both Lindblad and corotation resonances introduced through the 
vertical dependence. When E oc R~'^ and there is a central solar mass, they obtained: 



Results obtained from three dimensional calculations can be reproduced in two 
dimensional flat disk calculations if a gravitational softening parameter is adopted. 
When softening is introduced, the 1/R potential for the protoplanet is changed to 
1/ (i?^ + 6^)^''^ , the softening parameter being b. This can be regarded as accounting 
for the vertical averaging of the potential over the vertical thickness and accordingly b 
should be ~ H. 

We note that for their two dimensional calculations, Papaloizou & Larwood (2000) 
found that the torque results are dependent on the gravitational softening parameter 
used. Then tj and tg scale as 6^ '^^ and 6^'^, respectively, for b in the range O.AH-H. 
Papaloizou & Larwood (2000) find the approximate fits to tj and tg when there is a 
central solar mass and 7 = 3/2 to be: 




(76) 
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Figure 11. Density contour plots in a typical {x,y) plane for simulations with 
MpR^ / (Mi^H^) = 0.1 (top left panel), 0.3 (top right panel), 1.0 {bottom left panel) 
and 2.0 {bottom right panel). As MpR^ / {M-i,H^) increases, the wake becomes more 
prominent, material is accreted by the protoplanet and is pushed towards the radial 
boundaries as a gap is formed. The location of the characteristic ray defining the wake 
is also plotted (see Papaloizou et al 2004). The qualitative structure seen in these 
simulations is maintained once a quasi-steady state has been attained. (These plots 
are taken from Papaloizou et al 2004). 



Planet formation and migration 



54 




(77) 



and: 



te(yr) = 2.5 X 




(78) 

Here the gas disk has a mass Mqd contained within 5 AU and fg = 2.5b/H, where b is the 
gravitational softening length. Equations (76) and (77) agree when fg — 2.3, and 
give the characteristic timescale for migration of an Earth mass core in the minimum 
mass solar nebula as about 10^ yr. This timescale decreases linearly with mass, resulting 
in a survival problem for proto- Jovian mass planet cores. In this context, even the 
smallest cores considered of 5 have an infall timescale of only 2 x 10^ yr. 
The resolution of this problem may be provided by magnetic fields. When they are 
present and produce turbulence with density fluctuations that can be of order unity, 
type I migration may be stochastic (e.g.. Nelson & Papaloizou 2004, Laughlin et al 
2004). In this case, the migration direction of small mass cores becomes indeterminate 
over the timescales simulations can be performed and the final outcome may depend 
on long term global fluctuations in the MHD turbulence which in turn may depend 
on the global disk environment. We comment here that it is not clear that the disk 
behaves as a laminar disk with the time averaged density, upon which noisy fluctuations, 
with a characteristic time comparable to the orbital period and which obey Gaussian 
statistics, are superposed. If that were the case, it would be very difficult to affect 
the long term trend determined by the average laminar disk because of the large ratio 
of disk lifetime to orbital period. An additional recent investigation by Nelson (2005) 
suggests that turbulent fluctuations occur on all timescales up to the global diffusion 
or viscous time. Furthermore, the background perturbed by the planet is very different 
from laminar, with large density fluctuations encountering the planet and providing 
stochastic torques even when unperturbed by it. Consider a density fluctuation of 
order unity with length scale of order H a distance of order H from the planet. 
The characteristic speciflc torque acting on it is GT,R. Given the stochastic nature 
of turbulence, one might expect the specific torque acting on the planet to oscillate 
between ±GEit!. Note that this exceeds the net specific torque implied by equation (76) 
by a factor Nt ~ [{Mq/Mp){H/ R)^]{R/ H). Prom the discussion in section 4.2, the first 
factor should exceed unity for an embedded planet. Thus, such an object is inevitably 
subject to large torque fluctuations. The strength of the perturbation of the planet on 
the disk is measured by (Mp/ Mq){R/ H)^ . This perturbation only needs to produce a 
bias in the turbulent fluctuations resulting in a non zero mean, corresponding to the 
typical fluctuation reduced by a gactor Nf, for it to become comparable to the laminar 
type I value. 
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Also, global magnetic fields may be present and have significant effects on migration 
(Terquem 2003), or the disk may have long lived global distortions corresponding to non 
zero eccentricity (Papaloizou 2002) which may also significantly affect migration. The 
migration rates may also be affected by opacity transitions in the disk (Mcnou and 
Goodman 2004). Under some circumstances, significant fluctuations that can affect 
migration may occur in the non magnetic case for planets in the Saturn mass range 
(KoUer et al 2003) through the production of instabihties in inviscid gap edges. 

4-2.2. Type I migration in a turbulent disk: 

To illustrate stochastic migration when fully developed MHD turbulence is present 
in a disk with with H/R = 0.07 and no net magnetic flux, we show midplane density 
contours in figure 12, taken from Nelson & Papaloizou (2004). The presence of the 
protoplanet of mass 10 M® is just detectable, although the perturbation it makes to the 
disk is of lower amplitude than the perturbations generated by the turbulence. 
Midplane density contours for the same protoplanet in a laminar disk with identical 
parameters is given in figure 13. The disk response to the protoplanet is clearly much 
better defined in this case. 

The behavior of the torque per unit mass exerted by the disk on the protoplanet in 
the laminar case is illustrated in the left panel of figure 14, taken from Nelson & 
Papaloizou 2004. The upper line shows the torque due to the inner disk, the lower 
line shows the torque due to the outer disk, and the middle line shows the total torque. 
A well defined torque is produced, with an associated inward migration timescale in 
agreement with the standard linear theory. 

By way of contrast, the right panel of figure 14 shows the running time average of the 
torque per unit mass exerted by the disk in the turbulent case. The upper line is the 
running time average of the torque acting on the planet due to the inner disk. The 
lower line is that due to the outer disk. The middle (not straight) line is the running 
time average of the total torque. The straight line is the total torque exerted on the 
protoplanet in the laminar disk run. The total time averaged torque does not converge 
to a well defined value over the simulation run time of 20 orbits. 

These simulations indicate that, if the disk is not strictly laminar, type I migration 
may be significantly modified to the extent that the survival of protoplanetary cores 
may occur (this is also the case if the disk is magnetic or/and with gas not strictly 
on circular orbits). Furthermore, the enhanced mobility throughout the disk of a core 
undergoing stochastic migration may act to counteract the effects of core isolation (Rice 
& Armitage 2003). However, considerable theoretical developments need to be made in 
order to confidently assess the outcome of these phenomena. 
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Figure 12. This figure shows midplane density contours in the region close to the 
protoplanet. The protoplanet is located at {R,(j>)= (3,7r) in a disk with fully developed 
MHD turbulence with H/R = 0.07 and no net magnetic flux. The panels correspond 
to times 429.0687 {top left), 496.8698 {top right), 542.0710 {bottom left) and 634.4136 
{bottom right). The orbital period at i? = 1 is 27r in these units. Note that the presence 
of the protoplanet of mass 10 is just detectable, although the perturbation it makes 
to the disk is of lower amplitude than the perturbations generated by the turbulence. 
(These plots are taken from Nelson & Papaloizou 2004). 



Planet formation and migration 



57 



Laminar Disc - Mp=10 Eartti Masses 




2.5 - 



2.4 2.6 Z.8 3.0 3.2 3.4 3.6 3.8 



Figure 13. This figure shows midplane density contours for a laminar disk run 
with disk and protoplanet parameters otherwise identical to the simulation shown 
in figure 12. (This plot is taken from Nelson & Papaloizou 2004) 

4-3. Coorbital Torques and type III migration 

By coorbital torque, we mean a torque exerted on the protoplanet by disk material 
flowing through the orbit. This has been discussed recently by Masset & Papa- 
loizou (2003) and Artymowicz (2004) in a global context and by Papaloizou (2005) 
in a local shearing box context. 

One starts by considering a migrating protoplanet. When gap formation is only partial, 
material will flow through the orbit. As it does so, a torque will be exerted on the 
protoplanet. This may either act to assist the migration (positive feedback) or retard it 
(negative feedback). 

D'Angelo et al (2005) have recently studied coorbital torques acting on a Saturn mass 
planet in a massive disk under a similar set up to Masset & Papaloizou (2003). They 
used a multigrid system for which the resolution could be increased close to the planet, 
but not elsewhwere. They report only slow migration differing only marginally from 
type II migration. This was claimed to result from important torques acting within 
the Roche lobe. We firstly comment that their computational set up allowed mass 
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Figure 14. Left panel: Torque per unit mass exerted by tlie dislc on tlie protoplanet 
for tiie laminar disk simulation shown in figure 13. The upper line shows the torque 
due to the inner disk, the lower line shows the torque due to the outer disk, and the 
middle line shows the total torque. It is clear that a well defined torque is produced, 
with an associated migration time scale. 

Right panel: Running time average of the torque per unit mass exerted by the disk for 
the simulation shown in figure 12. The upper line is the running time average of the 
torque acting on the planet due to the inner disk. The lower line is that due to the outer 
disk. The middle (not straight) line is the running time average of the total torque. 
The straight line is the total torque exerted on the protoplanet in the comparable 
laminar disk run. We note that the total time averaged torque does not converge to a 
well defined value. (These plots are taken from Nelson & Papaloizou 2004). 

flow into the Roche lobe, effectively increasing the protoplanet mass by more than a 
factor of three. However, the increase in mass was not incorporated in the action of the 
protoplanet on the disk to determine its response, making the calculation inconsistent. 
Furthermore the results, checked with the highest resolution, shown in their figure 14, 
which differ from those used to draw final conclusions, indicate a fast inward migration 
time of ~ 600 orbits. This would not indicate important contributions from inside the 
Roche lobe as one would expect when there is not an increasing accumulation of mass 
there (see the appendix). However, the study of the migration of partially gap forming 
planets of modest mass, with a disk mass comparable to their own in the coorbital 
region, is difficult numerically and clearly further work on the problem is required. 
Masset & Papaloizou (2003) argue for positive feedback from coorbital torques because, 
for an outwardly migrating protoplanet, material traverses the coorbital zone from 
outside to inside. As it does so, it moves along the outer boundary of the coorbital 
zone occupied by material librating around the coorbital equilibrium (Lagrange) points 
passing close to the protoplanet. As this passage occurs, angular momentum is 
transferred to the protoplanet, a process which acts to assist the migration of the 
protoplanet and which accordingly gives a positive feedback. We now look at this 
process in more detail for an outwardly migrating protoplanet. Similar considerations 
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apply with appropriate sign reversals to an inwardly migrating protoplanet. 
The expected force exerted on the protoplanet by material flowing through the coorbital 
region is estimated as follows. The rate at which the mass flows through is 2'KRLdR/ dt. 
The outward momentum per unit mass imparted to the protoplanet when the disk 
matter moves across the coorbital region of width w is wuj/2, where uj is the rotation rate 
of the protoplanet. We adopt w — 2H, being appropriate for marginal gap formation. 
Then, the rate of transfer of momentum to the planet is Fcr — 27rRu;HT,(dR/dt). This 
can be expressed as: 

= 27tREuj^H^-^, (79) 

Cg Out' 

where Cs is the sound speed. Using ~ ujH, it may also be written in the form: 

F,, = (80) 

where — 'iwREH is the disk mass that would fill the coorbital zone of width 2H 
were it to do so at the background surface density. 

However, torques on the protoplanet do not only arise from material passing through the 
coorbital zone. Material that is forced to comove with the protoplanet, either because 
it has been accreted by it, or because it librates about coorbital equilibrium points, has 
to be acted on by the same force per unit mass as the protoplanet so as to maintain its 
migration. This results in an additional force acting on the protoplanet given by: 

^ dR 

F,rh=--MbUJ—, (81) 
where is the coorbital bound mass. The total force so far is thus: 

F,r = \{Ma-M^)uj^. (82) 

However, there are other forces. The forces acting due to density waves from the two 
sides will be affected by the migration and flow through and thus differ from the non 
migrating case. In this context, note that there is an associated asymmetry in the surface 
density profile. Thus, we should expect a wave torque component which is proportional 
to the migration speed. The indication from our simulations is that this acts as a drag 
on the protoplanet, as does M^. Accordingly, we shall consider the effects of asymmetric 
wave torques as modifying the value of Mi, and continue to use equation (82). 
Suppose now that the protoplanet is acted on by some external torque T^xt- The equation 
of motion governing the migration of the protoplanet of mass Mp with speed dR/dt, 
obtained by considering the conservation of angular momentum, is: 

1 dn 1 dn 

^MpRw^ = -{M,- M,)Ruj^ + Te,,. (83) 
Accordingly, we can consider the planet to move with an effective mass: 

Meff^Mp-{Md-Mb). (84) 
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The quantity (M^ — M^) has been called the coorbital mass deficit by Masset & 
Papaloizou (2003). When this is positive, there is an effective reduction in the inertia of 
the protoplanet. If there were no asymmetry in wave torques, the coorbital mass deficit 
would be the amount of mass evacuated from the gap region were it to be initially filled 
with the background density. Gap filling accordingly reduces the coorbital mass deficit. 
It is also clear that, because at least a partial gap is required, the protoplanet must be 
massive enough to produce a nonlinear response in the disk. 

Masset & Papaloizou (2003) and Papaloizou (2005) indeed found the coorbital mass 
deficit to be positive, resulting in positive feedback from coorbital torques. They also 
found that, in some circumstances, the coorbital mass deficit could become as large as 
the planet mass, so reducing the effective inertia to zero. Under these circumstances, 
a very fast or runaway migration may occur. Because the coorbital mass deficit is 
proportional to the disk surface density, the phenomenon requires a massive disk, 
typically an order of magnitude larger than the minimum mass solar nebula. 
An interesting feature of this type of migration is that it is most effective for protoplanets 
that bare marginally gap forming, and thus in the mass range 0.1-0.3 Mj. Such forming 
protoplanets may undergo a period of fast migration that brings them into the inner 
regions of the disk, where gap formation is easier, while they are still of sub- Jovian mass. 
It may thus form an explanation of why the 'hot Jupiters' are in general sub- Jovian 
in mass (see Masset & Papaloizou 2003 and Papaloizou & Nelson 2005). But note in 
this context that this explanation does not necessarily require type 111 migration, as 
migration rates tend to be maximal in the mass range 0.1-0.3 Mj for low mass disks too 
(see D'Angelo et al 2003) while lower masses may have migration inhibited by turbulence 
(Nelson & Papaloizou 2004). 

A graphical summary of the different types of migration is shown in figure 15. This 
shows the migration timescale as a function of protoplanet mass assuming a core of a 
constant mass of 15 M®. It allows for the possibihty of type I, type II and type III 
migration. It indicates that in order for a planet to form, type I migration, which 
produces radial infall timescales much shorter than the protoplanet growth time, has to 
be strongly suppressed relative to the prediction derived from a non magnetic, laminar 
and circular disk. However, beyond type 1 migration, survivability, at least to enable 
the formation of the planet, is not threatened on account of the shortened growth times. 
This is the case even if type III migration operates. 

5. Multiplanet systems and their interactions 

The presence of high orbital eccentricities amongst cxtrasolar planets is suggestive of 
a strong orbital relaxation or scattering process. For this to happen, formation must 
occur on a timescale short enough that strong dynamical interactions may take place 
subsequently. The gaseous environment of a disk may act to inhibit such interactions 
until it is removed. Gas free dynamical interactions of coplanar protoplanets formed on 
neighboring circular orbits have been considered by Weidenschilling & Marzari (1996) 
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Figure 15. This diagram shows the migration time estimated for a minimum mass 
solar nebula (solid line) and the mass growth time for protoplanet models with 15 
cores. The lowest runaway migration times, forming the lower boundary to the lower 
shaded region, are estimated for a model five times more massive than the minimum 
mass solar nebula. The dashed line shows the protoplanet growth time as a function of 
mass for models with standard opacity. The dashed-dotted line shows the growth time 
for models with 3 percent opacity. (This plot is taken from Papaloizou & Nelson 2005) 



and Rasio & Ford (1996). These may produce close scatterings and high eccentricities, 
but the observed distribution of extrasolar planets is not reproduced. 



5.1. Orbital relaxation in young planetary systems 

Papaloizou & Terquem (2001, see also Adams & Laughlin 2003) investigate a scenario 
in which 5 < < 100 planetary objects in the range of several Jupiter masses are 
assumed to form rapidly through fragmentation or gravitational instability occurring 
in a disk or protostellar envelope on a scale of Rmax = 100 AU. If these objects are 
put down in circular orbits about a solar mass star, at random in a volume contained 
within a spherical shell with inner and outer radii of O.lRmax and Rmax respectively, a 
strong relaxation on a timescale ~ 100 orbits occurs, leading to independence of details 
of initial conditions. 

For the range of values of A^ considered, the evolution is similar to that of a star 
cluster. Most objects escape, leaving at most 3 bound planets, the innermost with 
semi-major axis in the range 0. 01 Rmax^^-^ Rmax- However, close encounters or collisions 
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Figure 16. This figure shows the evolution of the semi-major axes {upper plot) and 
pericenter distances {lower plot) of the N = 8 planets in the system versus time 
(measured in units of Pq, being the orbital period at Rmax)- The lines correspond to 
the different planets. A line terminates just prior to the escape of a planet. (These 
plots are taken from Papaloizou & Terquem 2001) 

with the central star occur for about 10% of the cases. Tidal interaction, giving orbital 
circularization at fixed pericenter distance and leading to the formation of a very closely 
orbiting giant planet, is then a possibility. 

An example of a run with N = 8 masses selected uniformly at random in the interval 
(0, 5 X IO^^Mq) and with central stellar radius i?* = 1.337 x 10~^Rmax is illustrated in 
figure 16. At the end of this run, only 2 planets remain bound to the central star. 

The relaxation processes discussed above are more likely to apply to the more 
massive extrasolar planets exceeding 1 Mj. Their observed number increases with 
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distance (Zucker & Mazeh 2002 and section 2.2.1 above), as found in the simulations. 
So far, there are 23 isolated candidate objects with mass MpSini > 4.5 Mj, semi-major 
axes a in the range 0.3-3.9 AU and eccentricities e in the range 0.03-0.76. Amongst the 
'hot Jupiters' detected so far, r Boo is a possible candidate, being unusually massive 
with MpSini ~ 4 Mj. Note that, in multiple systems, only 3 planets have been detected 
with a mass larger than 4.5 Mj. They are actually significantly more massive than 
this hmit, being between 6 and 13 Mj. One of these systems contains a 7.7 Mj planet 
with a = 0.29 AU and e = 0.529, and a 16.9 Mj brown dwarf with a = 2.85 AU and 
e = 0.228. Such a system could well have formed through the process described above. 
The two other systems will be discussed below. 

5.2. Effect of an outer relaxing distribution of protoplanets on a low mass planet 
formed in an inner disc 

TcrqTicm & Papaloizou (2002) have considered the effects of an outer relaxing 
distribution of protoplanets on an inner disk in which one inner planet with a mass 
rUp = 0.3 Mj is formed on a timescale of 10^ yr. During the formation process, the 
eccentricity is assumed to be damped by tidal interaction with the disk while the planet 
is built up progressively with orbital radius a in the range 0.3-10 AU. 

An example of a run with N — 9 outer planets of mass 8 Mj and Rmax — 100 AU 
is illustrated in figure 17. After ~ 10^ yr, one outer planet with a ~ 8 AU and e ~ 0.6 
remains. The inner planet enters into a cycle in which e varies between 0.1 and 0.24. 
The mutual inchnation oscillates between and 30°. More extreme cycles have been 
produced in other examples. 

Thus, an outcome of an outer relaxing system could be an inner lower mass 
protoplanet with high orbital eccentricity. Among the observed planets, there are 
several candidates for which no companion has been detected so far. If the above 
scenario applies, there should be an outer massive planet with high eccentricity. 
Among the candidate systems with large radial velocity drifts potentially due to a 
companion selected by Fischer et al (2001) was HD 38529, which has a planet with 
mpSini = 0.76 Mj, a — 0.13 AU and e = 0.27. This is similar to the system illustrated 
in figure 17. More recently, the discovery of a companion with Mp sin i = 12.7 Mj, 
a = 3.7 AU and e = 0.36 has been announced (Fischer et al 2003). If sini = 0.8 is 
adopted for this system, a significant eccentricity is not excited in an initially circular 
inner planet orbit if the system is assumed coplanar. But, as illustrated in figure 18, 
for a mutual orbital inclination of 60°, a cycle is found in which the inner eccentricity 
reaches its observed value. 

As mentioned above, three of the observed multiple systems have a planet with a 
mass larger than 4.5 Mj. In two of these systems (HD 74156 and HD 38529), the massive 
planet has a semi-major axis on the order of 3.5 AU and a significant eccentricity. The 
companion has a lower mass, a semi-major axis in the range 0.1-0.3 AU and also a 
significant eccentricity. These systems may be explained by the scenario presented here. 
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Figure 17. Semi-major axes (in AU, upper plot) and eccentricities {lower plot) of the 
N = 9 outer planets {solid lines) and low-mass planet {dotted lines) versus time (in yr). 
Here, for the innermost low-mass planet planet, a = 0.3 AU initially. The eccentricity 
of the low-mass planet varies between 0.1 and 0.24 while a is almost constant. (These 
plots are taken from Terquem & Papaloizou 2002) 

6. Concluding renicirks 

The recent observations of extrasolar planets strongly suggest that gravitational 
interactions between forming planets and the nebula in which they are embedded form a 
key element in the understanding of the configurations of the systems. Given that the hot 
Jupiters have almost certainly not formed where they are observed today, they must have 
been subject to some orbital migration resulting from planet/disk tidal interaction. This 
type of interaction also leads to commensurable pairs of planets and apsidal resonances, 



Planet formation and migration 



65 



CD 
O 



o 



o 
o 




600 



Figure 18. Motion in HD 35829 for mutual inclinations of and 60°. The unchanging 
^ogio{a/ Rmax) for the two planets are given by the lower curves. The eccentricities 
by the upper curves. The oscillatory curve indicates a Kozai like cycle for the inner 
planet in the high inclination case. In the low inclination case the eccentricity of the 
inner planet remains close to zero. 



features that are also observed in some extrasolar systems. 

A substantial amount of work aimed at gaining a better understanding of 
planet/disk tidal interactions has been done in the last ten years. Three dimensional 
numerical simulations incorporating internally generated MHD turbulence rather than 
anomalous viscosity prescriptions have begun. They show that orbital migration is 
strongly affected by such processes, especially when the mass is in the Earth mass 
range. Analytical and numerical work has also shown that the presence of a large scale 
magnetic field or a finite orbital eccentricity can have a significant effect on migration. 

Strong gravitational interactions amongst a population of distant giant planets 
formed early in the life of the protostellar disk through gravitational fragmentation may 
also produce close orbiters and planets on highly eccentric orbits. The latter may belong 
to the original distribution or be formed by a later accumulation in an inner disk. 

We are currently at a stage where the important processes (gravitational instability, 
dust sedimentation, core accumulation, gas accretion and disk/planet interaction) can 
be studied in detail individually through numerical simulations. However, it is still not 
possible to take into account all of these processes acting simultaneously in a system with 
many embryos, to try to model the evolution of planetary systems from birth to final 
state. This would require very complex global simulations that cannot yet be performed. 
Also, for planetary systems, the evolutionary outcome, in the form of planetary masses 
and orbital configuration, depends very much on the initial conditions. The amount 
of mass present in the disk when planets begin to form is important in relation to 
gravitational instabilities, solid core masses that can be accumulated and migration 
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rates. Improved understanding will come as more powerful computers become available 
and observations put more constraints on realisable initial conditions. 
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APPENDIX 

Disk-protoplanet torques and angular momentum conservation 

We here consider a protoplanet in, either, a fixed circular orbit, or migrating in a 
circular orbit with slowly varying radius, that is embedded in and interacting with a 
gaseous disk. We show that under very general assumptions, the torque exerted between 
protoplanet and disk can be expressed as a combination of advective wave and viscous 
fluxes evaluated far from the protoplanet, together with a contribution accounting for 
the rate of change of angular momentum of the disk domain in the neighbourhod of the 
protoplanet when conditions are non stationary. In this way, the contributions relevant 
to different types of migration can be identified. 

In order to proceed, we need to consider the problem set up and distinguish between disk 
and protoplanet. It has been a common practice to treat the protoplanet as providing 
a forcing of the disk with a softened point mass potential. In that case, one may regard 
all the gas as belonging to the disk and identify the protoplanet mass as that associated 
with the potential, but only as long as the gas mass within a softening lengthscale 
of the protoplanet centre is negligible compared with that of the protoplanet. Another 
possibility would be to regard the protoplanet as residing interior to some domain within 
its Roche lobe. In that case, we consider the outer disk material to be unable to enter 
this domain or transmit angular momentum interior to it. Thus, the protoplanet itself 
does not change its mass or angular momentum through interacting with disk material 
at its boundary. But note that when this domain is well inside the Roche lobe, mass 
and angular momentum can still accumulate inside the Roche lobe but exterior to the 
protoplanet. 

Adopting cylindrical coordinates (r, ip, z) with origin at the central star of mass M^, we 
consider a protoplanet orbiting in a disk with uniform angular velocity r2k, where k is 
the unit vector in the vertical z direction, relative to an inertial frame. 
The basic equations governing the disk-protoplanet interaction are the conservation of 
mass and momentum in the form: 



dp 
dt 



V-(pv), 



(85) 
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and 

dv VP 

— + V • Vv + 21]kxv = f-V^o-V^i, (86) 

at p 

respectively. Here, p and v = {vr-iV^p, Vz) are the disk density and velocity, respectively. 
The pressure is P, V'o = —GM^j^r'^ + — Q?r'^ /2 is the gravitational potential due 
to the central object together with the centrifugal potential and -01 is the gravitational 
potential due to the orbiting protoplanet. 

The force per unit volume pf = p(/r, /<^, fz) is taken to be the divergence of a stress 
tensor, tij. Thus, in Cartesian coordinates r = (xi, X2, xs) for which f = {fi,f2,f3), 
pfi = dtij/dxj, where we have used the Einstein summation convention. Different effects 
may contribute to tij, including Navier Stokes viscosity in a laminar disk, self-gravity 
in a massive disk, or Maxwell stresses in a magnetised disk. 



Conservation of angular momentum: 

We may derive the relationship between the total torque acting between disk and 
protoplanet and wave and advected angular momentum fluxes in the disk by considering 
global angular momentum conservation. To do this, we start from the local conservation 
expressed by the azimuthal component of the equation of motion (86) in the form: 

On multiplying by pr and integrating over a flxed domain H, including the protoplanet 
but with a non interacting core cut out if needed, we obtain a statement of the total 
angular momentum conservation which, for an assumed steady state, takes the form: 

T = A{Fn, + Fstr). (88) 
Here, T is the total torque acting between protoplanet and disk: 

T = - / p^rdzdipdr. (89) 
Jh dip 

A{FRe + Fstr) denotes the difference of the total flux + Pstr (see below) evaluated 
on cylindrical bounding surfaces exterior to and interior to the protoplanet orbit, 
respectively. 

Pre is the anglar momentum flux associated with Reynolds stresses: 

Frs — J pVr{v^ + rQ)r^dzd(p, (90) 

with the integral being taken over a cylindrical bounding surface, S, assumed distant 
from the protoplanet. Fgtr is the contribution from the stress term: 

Fstr ^ ~ tr^r'^dzdip. (91) 
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The mass flow through a cyhndrical surface E', normal to a particular radius r in the 
disk, is given by: 

M = I pvrvdzdip. (92) 
Using this in equation (90), we obtain: 

FRe = pvr{v^ - v^y^dzdip + M{v^ + rVt)r. (93) 

Here the overline denotes azimuthal and vertical averaging with weight p. The first term 
can be identified as a wave or viscous flux (see for example Papaloizou and Lin 1984), as 
can Fsty. The second term may also be written as Mj, with j being the averaged specific 
angular momentum as viewed from the inertial frame. We shall assume the boundaries 
are far enough away from the protoplanet that j may be approximated as the Keplerian 
value. Note that, in a steady state, M is a constant. Further, turbulent or fluctuating 
states may be dealt with through appropriate time averaging and we assume this to be 
carried out where needed. Equation (88) then gives the torque exerted on the disk by 
the protoplanet as: 

T = MA(i)+A(F^e + i^str). (94) 

In this case, the torque can be expressed as the difference between boundary wave and 
viscous flux terms together with the difference between direct boundary output and 
input associated with the constant mass flux. Note that, when the protoplanet mass is 
reduced to zero, T = and we have MA(j) + A[FstrQ) = 0. Here Fgtro is the boundary 
viscous stress in the limit of vanishing protoplanet mass. This is just a statement 
that the advected angular momentum flow in a standard steady state disk is balanced 
by torques exerted through viscous, magnetic or gravitational stresses. Thus, when a 
protoplanet is present, we have: 

T^/^{Fne + Fstr-Fstro). (95) 

We see that, in this case, the torque exerted by the protoplanet on the disk can be 
expressed entirely in terms of boundary flux terms. 

A migrating protoplanet: 

To consider a migrating protoplanet, we suppose now that the angular velocity ^2 is a 
function of time, as is the orbital radius of the protoplanet R{t). We adopt a coordinate 
system such that the protoplanet appears stationary. Thus, this rotates with the time 
dependent angular velocity Q, while the other coordinate directions scale with the 
orbital radius R{t) . When conditions are quasi-steady, this corresponds to a self-similar 
solution. 

We define new dimensionles coordinates x[ — Xi/R{t). We also adopt a scaled 
dimensionless time coordinate r, which is such that dr — fldt. 
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Further, we adopt new dimensionless velocity components given by: 




V r dR 



(96) 



m E?n dt ■ 

In terms of the new variables, the continuity equation can be written: 




(97) 



where ps = pR^/Rl, where Rq is a constant with the dimensions of length that can be 
taken to be the value of R at some time t — to that can be chosen for convenience and 
the prime in V' denotes differentiation with respect to x^. 
Similarly, the eqation of motion becomes: 



Note that we here assume that Q is a function of R. In fact, we shall use the Keplerian 
form fl^ = GMi,/R^. In addition, we shall adopt dR/dt = eRQ, where e is a constant 
of small magnitude corresponding to a slow radial migration. In that case, when the 
pressure and f can be scaled appropriately, it is possible to seek a steady solution of 
equations (97)-(98) for the primed state variables that is independent of r. When the 
equation of state is P = pc^, with being the sound speed, assumed to be oc r~^/^, 
the stress tensor corresponding to f scales as P, and the gravitational softening length 
scales as it! if it is being used, R may be scaled out of the basic equations thus enabling 
a steady state to be sought. 

We comment that the last two terms on the left hand side of equation (98) arise from the 
transformation to an expanding/contracting rotating frame. In this frame, they appear 
as the additional forces/torques required to induce the apparent modified disk flow. 

Conservation of angular momentum: 

As in the fixed orbit case, we formulate the global conservation of angular momentum 
derived from equation (98). On multiplying by p^r' and integrating over a fixed 
domain V, which is similar to Ti but is defined in the primed coordinates, the azimthal 
component of equation (98) gives the total angular momentum conservation in the form: 



i?2Q Dt 



1 D 



(pW) +2(kxv') + 



(kxr)/ djR^Q) dR r' dPR 
Q2p2 dR Itt^ Q?R df^ 




(98) 




(99) 




(100) 



Planet formation and migration 



70 



The angular momentum content within T>, as seen in the inertial frame, is: 

J = R^nii^ ps{v'^ + r'y^dr'dz'dip', (101) 
while the anglar momentum flux associated with Reynolds stresses is: 

F«e = R'^^'^Rl I Psv'rK + rVdz'dip', (102) 

with the integral being taken over a bounding surface S, assumed distant from the 
protoplanet amd fixed in the primed coordinates. The contribution from the stress 
term is: 

F,tr = -R^ [ tr' u,'r'^dz'dip'. (103) 
Js 

As before, the mass flow through a cylindrical surface normal to a particular radius, r', 
in the disk which remains fixed in the primed coordinates, is given by: 

M^nRl ( p^v'j'dz'd^'. (104) 
Js 

Using this in equation (102), we obtain: 

FRe = R^^'^Rl [ P3v',{v'^ - v^y^dz'dif' + MR^n{vJ + r')r'. (105) 

J s 

Here the overline denotes azimuthal and vertical averaging with weight p^. The first term 
can be identified as a wave or viscous flux, as can Fgtr- The second term may also be 
written as Mj, with j being the averaged speciflc angular momentum. Again we assume 
the boundaries to be sufficiently far from the protoplanet that j may be approximated 
as Keplerian. In general, M is a function of both r' and r or time. However, in a 
putative steady state, as seen in the corotating comoving coordinate system (primed 
variables independent of r), it is a function of time only. 

Using the above, equation (99) gives the torque exerted on the disk by the protoplanet 
as: 

T = + A(Mj) + AF^„ (106) 
or 

where we have combined the viscous and wave fiuxes in Thus the torque can be 
expressed as the difference between boundary wave and viscous flux terms together with 
the difference between direct output and input associated with the mass throughput in 
combination with the rate of growth of the angular momentum contained in the domain 
T>. It is the latter two terms that are responsible for the coorbital torques. When the 
angular momentum change associated with mass throughput is not accounted for by an 
increase or decrease in the angular momentum content of the domain traversed, it is 
accounted for by the protoplanet disk torque, which means that there has to be angular 
momentum exchange with the protoplanet. 

Let us look at this more closely. Although conclusions do not depend on this, for 
simplicity we consider an assumed steady state self-similar case, then M does not vary 
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in space. But note that, from equation (101), because of the dependence on R and Q, 
VtdJ/dT = J/{2R){dR/dt) ^ 0. In this case we have: 

T ^ —— + M{Aj) + AF^^. (107) 

Suppose now viscous effects are small, so that M as viewed from the comoving coordinate 
frame we use is due to the flow of matter induced by the migration. When the mass 
of the protoplanet is negligible, but the migration maintained, we expect that, as then 
there are no excited waves and T — 0: 

^f + M(A,) = 0. (108) 

Here Jq is the angular momentum content of T> in this limit. Using this in equation (108) 
gives: 

^ (J- Jo)dR , ^ , 
Now, conservation of angular momentum of the planetary orbit gives: 

r = (110) 

2Rdt' ^ ^ 

where Jp is the angular momentum of the protoplanet. Thus we have: 

a-(Jo-J))jg^ 

2R dt ^ ' 

Thus, as far as the orbital evolution is concerned, the protoplanet orbital angular 
momentum is reduced by the coorbital angular momentum deficit Jq — J. Other things 
being equal, this is expected to be positive if the protoplanet generates a partial gap, 
and in that case, once the effect of any asymmetric wave torques is taken into account, 
is clearly related to the coorbital mass deficit considered in section (4.3. There, it 
leads to a positive feedback enhancing the migration. On the other hand, if Jq — J is 
negative on account of excess material accreted in the vicinity of the planet increasing 
the angular momentum content interior to the Roche lobe, then the reverse effect of 
negative feedback occurs. 

Note further that, if the structure interior to the Roche lobe is quasi steady, and if the 
total mass of material there, not counted as part of the protoplanet, is small compared 
to that of the protoplanet, only a small contribution to the disk-protoplanet torque is 
expected from these regions. 
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